En esta práctica se va a analizar un dataset que aborda el rendimiento en la asignatura Matemáticas de los alumnos en educación secundaria de dos centros escolares portugueses. Los datos se han recopilado mediante el uso de informes y cuestionarios escolares. Dicho dataset se puede descargar mediante el siguiente enlace:

https://archive.ics.uci.edu/dataset/320/student+performance

En primer lugar, se realizará un análisis exploratorio previo de los datos para identificar posibles valores perdidos y valores extremos o “outliers”, y se tomarán las decisiones correspondientes para tratarlos.

En segundo lugar, se realizará un Análisis de Componentes Principales (ACP), esto es, se va a condensar la información aportada por las variables originales en unas pocas combinaciones lineales de ellas, con el objetivo de hacer una reducción de la dimensión. Estas se buscan con máxima varianza y perpendiculares entre sí, de forma que cada una sigue la dirección en la que las observaciones varían más y no está correlacionada con las anteriores.

En tercer lugar, se realizará un Análisis Factorial (AF), esto es, se van a identificar variables latentes (no observables) con una alta correlación con un grupo de variables observables y correlación prácticamente nula con el resto. Así, se hará una reducción de la dimensión.

En cuarto lugar, se realizará un Análisis Discriminante Lineal (ADL) y otro Cuadrático (ADC), verificando las hipótesis necesarias de normalidad. El Análisis Discriminante es un método de clasificación de variables cualitativas que permitirá clasificar nuevas observaciones según sus características (variables explicativas o predictores) en las distintas categorías de la variable cualitativa respuesta.

Finalmente, se realizará un Análisis Cluster (AC). El Análisis Cluster es una técnica multivariante cuyo principal objetivo es agrupar objetos formando conglomerados (clusters) con un alto grado de homogeneidad interna y heterogeneidad externa. En otras palabras, el AC es un procedimiento exploratorio que permite encontrar estructuras con similitudes en un conjunto de datos con cierta variabilidad.

1 Cargando Paquetes y Leyendo el Conjunto de Datos

1.1 Cargar e instalar paquetes de R

El siguiente módulo de código fuente se encarga de cargar, si ya están instalados, todos los paquetes que se utilizarán en esta sesión de R. Si bien un paquete R se puede cargar en cualquier momento cuando se vaya a utilizar, es recomendable optimizar sus llamadas con este fragmento de código al principio.

Cargar un paquete en una sesión de R requiere que ya esté instalado. Si no es así, el primer paso es ejecutar la sentencia:

install.packages(“name_of_the_library”)

#########################################
# Loading necessary packages and reason #
#########################################

# This is an example of the first installation of a package
# Only runs once if the package is not installed
# Once it is installed this sentence has to be commented (not to run again)
# install.packages("summarytools")

# Package required to call 'read_xlsx' function (loading '.xlsx' data format)
library(readxl)

# Package required to call 'freq' and 'descr' functions (descriptive statistics)
library(summarytools)

# Package required to call 'factanal' function
library(stats)

# Package required to call 'stat.desc' function
library(pastecs)

# Package required to call 'cortest.bartlett' function
library(psych)

# Package required to call 'hetcor' function
library(polycor)

# Package required to call 'ggcorrplot' function
library(ggcorrplot)

# Package required to call 'corrplot' function
library(corrplot)

# Package required to call 'rplot' function
library(corrr)

# Package required to call 'ggplot' function (graphical tools)
library(ggplot2)

# Package required to call 'ggarrange' function (graphical tools)
library(ggpubr)

# Package required to call 'fviz_pca_var, fviz_pca_ind and fviz_pca' functions
library(factoextra)

# Package required to call 'scatterplot3d' function
library(scatterplot3d)

# Package required to call 'mutate' function
library(tidyverse)

# Package required to call 'clusGap' function
library(cluster)

# Package required to call 'ggdendrogram' function
library(ggdendro)

# Package required to call 'grid.arrange' function
library(gridExtra)

# Package required to call 'melt' function
library(reshape2)

# Package required to call 'mvn' function
library(MVN)

# Package required to call 'boxM' function
library(biotools)

# Package required to call 'partimat' function
library(klaR)

# Package required to call 'summarise' function
library(dplyr)

# Package required to call 'createDataPartition' function
library(caret)

1.2 Descripción del conjunto de datos student_mat_DB

El fichero student_mat_DB contiene datos recogidos de 395 estudiantes recogidos en 31 variables. A continuación se muestra una descripción de cada variable del conjunto de datos:

  • Centro escolar del estudiante (school)
  • Sexo del estudiante (sex)
  • Edad del estudiante (age)
  • Tipo de dirección de vivienda del estudiante (address)
  • Número de miembros de la familia del estudiante (famsize)
  • Estado de convivencia de los padres del estudiante (Pstatus)
  • Educación de la madre del estudiante (Medu)
  • Educación del padre del estudiante (Fedu)
  • Trabajo de la madre del estudiante (Mjob)
  • Trabajo del padre del estudiante (Fjob)
  • Motivo del estudiante por el cual eligió su centro escolar (reason)
  • Tutor del estudiante (guardian)
  • Tiempo de viaje de casa del estudiante al centro escolar (traveltime)
  • Tiempo de estudio semanal del estudiante (studytime)
  • Número de materias suspendidas en el pasado (failures)
  • El estudiante recibe apoyo educativo adicional (schoolsup)
  • El estudiante recibe apoyo educativo familiar (famsup)
  • El estudiante recibe clases extras pagadas de la asignatura Matemáticas (paid)
  • El estudiante está apuntado en actividades extracurriculares (activities)
  • El estudiante asistió de pequeño a la guardería (nursery)
  • El estudiante quiere cursar estudios superiores (higher)
  • El estudiante tiene acceso a Internet en casa (internet)
  • El estudiante se encuentra en una relación de pareja (romantic)
  • Calidad de las relaciones familiares del estudiante (famrel)
  • Cantidad de tiempo libre que posee el estudiante después del colegio (freetime)
  • Frecuencia con la que el estudiante sale a la calle con sus amigos (goout)
  • Cantidad de consumo de alcohol del estudiante en jornada laboral (Dalc)
  • Cantidad de consumo de alcohol del estudiante en el fin de semana (Walc)
  • Estado de salud actual del estudiante (health)
  • Número de veces que el estudiante ha faltado al colegio (absences)
  • Nota final del estudiante en la asignatura Matemáticas (G3)

1.2.1 Lectura del conjunto de datos

El archivo de datos tiene extensión .xlsx (Microsoft Excel). Para cargar el fichero se utiliza la función read_excel() dentro del paquete readxl.

Para cargar la base de datos, la sesión de R debe ejecutarse en el mismo directorio en el que se encuentra el fichero. Para ello, se puede usar la función setwd(“Ubicación del archivo de datos”) o se puede pulsar Session/Set Working Directory/Choose Directory en la barra de menús de RStudio. Asumiendo que la sesión de R está correctamente direccionada, el siguiente código carga el fichero de datos en un data.frame.

# Setting the work directory
setwd("C:/Users/LENOVO/Desktop/EMV/practicas/practica_final")

# Loading a .xlsx (excel) file.
# The output of this function is already a data.frame object
# Remember that package 'readxl' is required
data_xlsx<-read_excel("student_mat_DB.xlsx", sheet = 2)

# This sentence identifies the type of object that the identifier represents
class(data_xlsx)
## [1] "tbl_df"     "tbl"        "data.frame"

2 Análisis Explotario Univariante

2.1 Recodificaciones o agrupaciones de datos

Se ha recodificado el conjunto de datos original. Esta recodificación ya se encuentra realizada y detallada en el archivo student_mat_DB.xlsx.

Cabe destacar que es cierto que, por el carácter de las variables, principalmente categóricas pero recodificadas, podrían discutirse los resultados obtenidos tras aplicar las técnicas para la reducción de la dimensión (ACP y AF, los cuales realizaremos posteriormente) puesto que estos trabajan con variables cuantitativas, pero también es cierto que si las variables son binarias (como ocurre en varias variables de nuestro conjunto de datos), todavía se puede justificar su uso y usabilidad.

Hay que tener cuidado con las variables codificadas que eran categóricas originalmente y que tienen más de dos niveles de respuesta porque quizá no interese incluirlas. Para decidir si las incluímos o no, tendría que ser el que ha diseñado el dataset el que nos dijera si esas variables son relevantes para el problema bajo estudio. Asumiendo que lo son, una forma de ver si pueden ser útiles es hacer un análisis descriptivo de ellas y ver si los niveles están compensados (más o menos hay el mismo número de individuos en cada uno). Si lo están es una variable que podría aportar información, pero si hay algún nivel muy bajo o si todo se acumula en un nivel, habría que plantearse descartar la variable. No hay una forma definitiva de actuar, depende de lo que busquemos, pero este camino podría ayudar.

2.1.1 Variable Medu

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Medu))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0      3      0.76           0.76      0.76           0.76
##           1     59     14.94          15.70     14.94          15.70
##           2    103     26.08          41.77     26.08          41.77
##           3     99     25.06          66.84     25.06          66.84
##           4    131     33.16         100.00     33.16         100.00
##        <NA>      0                               0.00         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p1<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Medu)))+geom_bar()+
  coord_polar("y")+labs(x="Medu", y="%")

p2<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Medu)))+geom_bar()+
  labs(x="Medu", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p1, p2, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que la frecuencia del nivel asociado al valor 0 es muy baja (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.

# The variable 'Medu' is eliminated
data_xlsx<-data_xlsx[, -c(7)]

2.1.2 Variable Fedu

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Fedu))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0      1      0.26           0.26      0.25           0.25
##           1     77     20.26          20.53     19.49          19.75
##           2    114     30.00          50.53     28.86          48.61
##           3     97     25.53          76.05     24.56          73.16
##           4     91     23.95         100.00     23.04          96.20
##        <NA>     15                               3.80         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p3<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fedu)))+geom_bar()+
  coord_polar("y")+labs(x="Fedu", y="%")

p4<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fedu)))+geom_bar()+
  labs(x="Fedu", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p3, p4, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que la frecuencia del nivel asociado al valor 0 es muy baja (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.

# The variable 'Fedu' is eliminated
data_xlsx<-data_xlsx[, -c(7)]

2.1.3 Variable Mjob

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Mjob))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0     55     14.29          14.29     13.92          13.92
##           1     33      8.57          22.86      8.35          22.28
##           2    100     25.97          48.83     25.32          47.59
##           3     59     15.32          64.16     14.94          62.53
##           4    138     35.84         100.00     34.94          97.47
##        <NA>     10                               2.53         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p5<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Mjob)))+geom_bar()+
  coord_polar("y")+labs(x="Mjob", y="%")

p6<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Mjob)))+geom_bar()+
  labs(x="Mjob", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p5, p6, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que más o menos están compensados los niveles de esta variable. Por tanto, mantenemos dicha variable en nuestro conjunto de datos.

2.1.4 Variable Fjob

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Fjob))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0     29      7.34           7.34      7.34           7.34
##           1     18      4.56          11.90      4.56          11.90
##           2    111     28.10          40.00     28.10          40.00
##           3     20      5.06          45.06      5.06          45.06
##           4    217     54.94         100.00     54.94         100.00
##        <NA>      0                               0.00         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p7<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fjob)))+geom_bar()+
  coord_polar("y")+labs(x="Fjob", y="%")

p8<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fjob)))+geom_bar()+
  labs(x="Fjob", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p7, p8, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que la frecuencia del nivel asociado al valor 4 es muy alta (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.

# The variable 'Fjob' is eliminated
data_xlsx<-data_xlsx[, -c(8)]

2.1.5 Variable reason

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$reason))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0    109     27.59          27.59     27.59          27.59
##           1    105     26.58          54.18     26.58          54.18
##           2    145     36.71          90.89     36.71          90.89
##           3     36      9.11         100.00      9.11         100.00
##        <NA>      0                               0.00         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p9<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$reason)))+geom_bar()+
  coord_polar("y")+labs(x="reason", y="%")

p10<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$reason)))+geom_bar()+
  labs(x="reason", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p9, p10, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que más o menos están compensados los niveles de esta variable. Por tanto, mantenemos dicha variable en nuestro conjunto de datos.

2.1.6 Variable guardian

# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$guardian))
## Frequencies  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0    273     69.11          69.11     69.11          69.11
##           1     90     22.78          91.90     22.78          91.90
##           2     32      8.10         100.00      8.10         100.00
##        <NA>      0                               0.00         100.00
##       Total    395    100.00         100.00    100.00         100.00
# Pie chart and bar graph
p11<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$guardian)))+geom_bar()+
  coord_polar("y")+labs(x="guardian", y="%")

p12<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$guardian)))+geom_bar()+
  labs(x="guardian", y="%")

# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p11, p12, nrow = 1, ncol=2, common.legend = TRUE)

Vemos que la frecuencia del nivel asociado al valor 0 es muy alta (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.

# The variable 'guardian' is eliminated
data_xlsx<-data_xlsx[, -c(9)]

2.2 Explorando el conjunto de datos

# Retrieving the name of all variables
colnames(data_xlsx)
##  [1] "school"     "sex"        "age"        "address"    "famsize"   
##  [6] "Pstatus"    "Mjob"       "reason"     "traveltime" "studytime" 
## [11] "failures"   "schoolsup"  "famsup"     "paid"       "activities"
## [16] "nursery"    "higher"     "internet"   "romantic"   "famrel"    
## [21] "freetime"   "goout"      "Dalc"       "Walc"       "health"    
## [26] "absences"   "G3"
# Displaying a few records
head(data_xlsx, n=10)
data_xlsx
# Checking the type of the variables
sapply(data_xlsx, class)
##     school        sex        age    address    famsize    Pstatus       Mjob 
##  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric" 
##     reason traveltime  studytime   failures  schoolsup     famsup       paid 
##  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric" 
## activities    nursery     higher   internet   romantic     famrel   freetime 
##  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric" 
##      goout       Dalc       Walc     health   absences         G3 
##  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"  "numeric"

2.3 Valores perdidos (NA)

# Construction of the function that calculates the percentage of missing values for each variable
porcentaje_NA<-function(data){
  c=(sum(is.na(data)))/length(data)*100
  return(c)
}

# Checking for missing values
cont_NA<-apply(data_xlsx, 2, porcentaje_NA)
cont_NA
##     school        sex        age    address    famsize    Pstatus       Mjob 
##   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   2.531646 
##     reason traveltime  studytime   failures  schoolsup     famsup       paid 
##   0.000000   0.000000  13.670886   0.000000   0.000000   0.000000   0.000000 
## activities    nursery     higher   internet   romantic     famrel   freetime 
##   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000 
##      goout       Dalc       Walc     health   absences         G3 
##   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000

Se observa que la variable Mjob tiene menos de un 5% de valores perdidos, a diferencia de la variable studytime, que sí que tiene más de un 5%. Por lo tanto, hay que analizar en dicha variable studytime si el patrón seguido por dichos valores perdidos es aleatorio o no. Para ello, se estudia la homogeneidad según grupos (NA y no NA) con otras variables según un Test de Student (ya que la variable es continua).

Se define una nueva variable studytime_NA que toma el valor 0 si el correspondiente dato no es perdido y el valor 1 si el correspondiente dato es perdido. A continuación, se realiza un test de contraste de igualdad de medias entre los grupos de la variable traveltime (que no tiene valores perdidos) correspondientes al 0 y al 1 de estas nuevas variables.

# Construction of the function that replaces missing values of a variable with 1 and non-missing values with 0
indicadora_NA<-function(data, na.rm=F){
  data[!is.na(data)]<-0
  data[is.na(data)]<-1
  return(data)
}

studytime_NA<-indicadora_NA(data_xlsx$studytime)

# Test de Student
t.test(data_xlsx$traveltime~studytime_NA, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  data_xlsx$traveltime by studytime_NA
## t = -1.2192, df = 393, p-value = 0.2235
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -0.32519085  0.07624983
## sample estimates:
## mean in group 0 mean in group 1 
##        1.431085        1.555556

Como el p-valor es grande en todos los casos (p-valor > 0.15), no se encuentran evidencias para rechazar la hipótesis nula de igualdad de medias, de forma que se acepta la hipótesis de homogeneidad y se concluye que el patrón es aleatorio. En este caso, como las variables son cuantitativas, la decisión para los datos perdidos es reemplazarlos por la media de su variable.

# Construction of the function that handles missing values
not_available<-function(data, na.rm=F){
  data[is.na(data)]<-mean(data, na.rm=T)
  return(data)
}

data_xlsx$Mjob<-not_available(data_xlsx$Mjob)
data_xlsx$studytime<-not_available(data_xlsx$studytime)

# Viewing the data again
head(data_xlsx, n=10)
# Checking for missing values
no_cont_NA<-apply(data_xlsx, 2, porcentaje_NA)
no_cont_NA
##     school        sex        age    address    famsize    Pstatus       Mjob 
##          0          0          0          0          0          0          0 
##     reason traveltime  studytime   failures  schoolsup     famsup       paid 
##          0          0          0          0          0          0          0 
## activities    nursery     higher   internet   romantic     famrel   freetime 
##          0          0          0          0          0          0          0 
##      goout       Dalc       Walc     health   absences         G3 
##          0          0          0          0          0          0

2.4 Análisis descriptivo numérico clásico

A continuación, se realiza un análisis descriptivo de los datos. La función stat.desc() del paquete pastecs proporciona las medidas de posición, dispersión y forma más importantes:

  • Número de valores (nbr.val)
  • Número de valores nulos (nbr.null)
  • Número de valores perdidos (nbr.na)
  • Valor mínimo de la variable (min)
  • Valor máximo de la variable (max)
  • Rango de valores de la variables (max-min) (range)
  • Suma de los valores no perdidos (sum).
  • Mediana (median).
  • Media (mean).
  • Error estándar de la media (SE.mean)
  • Intervalo de confianza de la media (CI.mean)
  • Varianza (var)
  • Desviación típiaca (std.dev)
  • Coeficiente de variación (coef.var)
  • Coeficiente de asimetría (skewness)
  • Estadístico para contrastar si el coeficiente de asimetría es cero (skew.2SE)
  • Coeficiente de curtosis (kurtosis)
  • Estadístico para contrastar si el coeficiente de curtosis es cero (kurt.2SE).
  • Estadístico de Shapiro-Wilks (normtest.W)
  • Probabilidad asociada al estadístico de Shapiro-Wilks (normtest.p)
# Basic descriptive statistics
stat.desc(data_xlsx, norm=TRUE)

Para más detalles sobre esta función se puede consultar el siguiente enlace:

https://www.rdocumentation.org/packages/pastecs/versions/1.3.21/topics/stat.desc

Como esta función no muestra los cuartiles, se utiliza la función summary() para obtenerlos junto al resto de medidas de posición:

# Basic descriptive statistics
summary(data_xlsx)
##      school            sex              age          address      
##  Min.   :0.0000   Min.   :0.0000   Min.   :15.0   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:16.0   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :17.0   Median :0.0000  
##  Mean   :0.1165   Mean   :0.4734   Mean   :16.7   Mean   :0.2228  
##  3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:18.0   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :22.0   Max.   :1.0000  
##     famsize          Pstatus            Mjob           reason     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:0.000  
##  Median :1.0000   Median :0.0000   Median :2.499   Median :1.000  
##  Mean   :0.7114   Mean   :0.1038   Mean   :2.499   Mean   :1.273  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:4.000   3rd Qu.:2.000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :4.000   Max.   :3.000  
##    traveltime      studytime        failures        schoolsup     
##  Min.   :1.000   Min.   :1.000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:1.000   1st Qu.:2.000   1st Qu.:0.0000   1st Qu.:1.0000  
##  Median :1.000   Median :2.000   Median :0.0000   Median :1.0000  
##  Mean   :1.448   Mean   :2.023   Mean   :0.3342   Mean   :0.8709  
##  3rd Qu.:2.000   3rd Qu.:2.023   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :4.000   Max.   :4.000   Max.   :3.0000   Max.   :1.0000  
##      famsup            paid          activities        nursery      
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.0000   Median :1.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.3873   Mean   :0.5418   Mean   :0.4911   Mean   :0.2051  
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##      higher           internet         romantic          famrel     
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :1.000  
##  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:4.000  
##  Median :0.00000   Median :0.0000   Median :1.0000   Median :4.000  
##  Mean   :0.05063   Mean   :0.1671   Mean   :0.6658   Mean   :3.944  
##  3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:5.000  
##  Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :5.000  
##     freetime         goout            Dalc            Walc      
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:3.000   1st Qu.:2.000   1st Qu.:1.000   1st Qu.:1.000  
##  Median :3.000   Median :3.000   Median :1.000   Median :2.000  
##  Mean   :3.235   Mean   :3.109   Mean   :1.481   Mean   :2.291  
##  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:2.000   3rd Qu.:3.000  
##  Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.000  
##      health         absences            G3       
##  Min.   :1.000   Min.   : 0.000   Min.   : 0.00  
##  1st Qu.:3.000   1st Qu.: 0.000   1st Qu.: 8.00  
##  Median :4.000   Median : 4.000   Median :11.00  
##  Mean   :3.554   Mean   : 5.709   Mean   :10.42  
##  3rd Qu.:5.000   3rd Qu.: 8.000   3rd Qu.:14.00  
##  Max.   :5.000   Max.   :75.000   Max.   :20.00

2.5 Outliers

Para detectar los valores atípicos (outliers), hay que hacer un análisis exploratorio gráfico previo, construyendo los boxplots de todas las variables.

# Boxplots of all variables together
# This visualization is not the best due to the difference between the scales
boxplot(data_xlsx, main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(data_xlsx)))

Se ha tomado la decisión de sustituir los outliers por la media de sus variables. Para ello, se ha construido la función outlier(), que realiza la detección y manipulación de los mismos. Cabe destacar que la decisión de sustituir los outliers por la media requeriría un análisis previo de la causa de estos valores atípicos.

# Recursive function that modifies outliers by the mean of their variable
outlier<-function(data, na.rm=T){
  H<-1.5*IQR(data)
  data[data<quantile(data, 0.25, na.rm = T)-H]<-NA
  data[data>quantile(data, 0.75, na.rm = T)+H]<-NA
  data[is.na(data)]<-mean(data, na.rm = T)

  H<-1.5*IQR(data)

  if (TRUE %in% (data<quantile(data, 0.25, na.rm = T)-H) |
      TRUE %in% (data>quantile(data, 0.75, na.rm = T)+H))
    outlier(data)
  else
    return(data)
}

# This data.frame is to preserve original data once the outliers are modified
original_data_xlsx<-data_xlsx

# Called to outlier function for each variable identified with outliers
data_xlsx<-as.data.frame(apply(data_xlsx, 2, outlier))

# We compare the original data and the fixed ones with respective boxplots 
par(mfrow=c(1,2))
# Boxplot original data
boxplot(original_data_xlsx, main="Original data",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(original_data_xlsx)))
# Boxplot fixed data
boxplot(data_xlsx, main="Data with no outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(data_xlsx)))

Vemos que se aprecian muy bien los boxplots debido a los altos valores de la penúltima variable (absences) en comparación con el resto. A continuación se repite el análisis de outliers pero con los datos normalizados para que se pueda tener una mejor visualización de las gráficas anteriores.

# Normalized original data
normalized_data_xlsx<-original_data_xlsx
normalized_data_xlsx<-scale(original_data_xlsx)

# Boxplots of all variables together
# This visualization is not the best due to the difference between the scales
boxplot(normalized_data_xlsx, main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(normalized_data_xlsx)))

# This data.frame is to preserve original data once the outliers are modified
normalized_original_data_xlsx<-normalized_data_xlsx

# Called to outlier function for each variable identified with outliers
normalized_data_xlsx<-as.data.frame(apply(normalized_data_xlsx, 2, outlier))

# We compare the normalized original data and the fixed ones with respective boxplots 
par(mfrow=c(1,2))
# Boxplot normalized original data
boxplot(normalized_original_data_xlsx, main="Normalized original data",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(normalized_original_data_xlsx)))
# Boxplot fixed data
boxplot(normalized_data_xlsx, main="Normalized data with no outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=c(1:ncol(normalized_data_xlsx)))

3 Análisis de Componentes Principales (ACP)

3.1 Correlación

En primer lugar es necesario comprobar que las variables no son independientes. A nivel de la muestra recogida en la base de datos, esto se puede hacer calculando y observando la matriz de correlaciones. A nivel poblacional se justificará que hay correlación realizando el test de Bartlett (el contraste de esfericidad de Bartlett permite comprobar si las correlaciones son distintas de 0 de modo significativo. La hipótesis nula es que \(det(R)=1\) , siendo \(R\) la matriz de correlaciones). El siguiente código realiza las dos comprobaciones indicadas.

###############################
# Correlation at sample level #
###############################

# Are the variables correlated at sample level?
correlation_matrix<-cor(data_xlsx)
correlation_matrix
##            school          sex         age address       famsize Pstatus
## school          1           NA          NA      NA            NA      NA
## sex            NA  1.000000000 -0.04064933      NA -0.0898618018      NA
## age            NA -0.040649328  1.00000000      NA -0.0455885843      NA
## address        NA           NA          NA       1            NA      NA
## famsize        NA -0.089861802 -0.04558858      NA  1.0000000000      NA
## Pstatus        NA           NA          NA      NA            NA       1
## Mjob           NA -0.102773909  0.08962071      NA  0.0915722328      NA
## reason         NA  0.009821817 -0.01340060      NA  0.0009812261      NA
## traveltime     NA  0.020908545  0.11516528      NA -0.0577081482      NA
## studytime      NA  0.187170093 -0.05072535      NA  0.0341291336      NA
## failures       NA           NA          NA      NA            NA      NA
## schoolsup      NA           NA          NA      NA            NA      NA
## famsup         NA  0.151622617  0.13021861      NA -0.1128930671      NA
## paid           NA  0.129125619  0.02681464      NA -0.0138821617      NA
## activities     NA -0.099833177  0.09440942      NA -0.0001131763      NA
## nursery        NA           NA          NA      NA            NA      NA
## higher         NA           NA          NA      NA            NA      NA
## internet       NA           NA          NA      NA            NA      NA
## romantic       NA  0.102023009 -0.15316341      NA  0.0343948162      NA
## famrel         NA  0.069375471  0.03237037      NA  0.0069223355      NA
## freetime       NA  0.207351483  0.02663625      NA -0.0131216979      NA
## goout          NA  0.075897396  0.11147582      NA -0.0230643843      NA
## Dalc           NA  0.224996943  0.09032711      NA -0.1485945736      NA
## Walc           NA  0.274193767  0.09719130      NA -0.1034250076      NA
## health         NA  0.143588173 -0.04372758      NA  0.0289916591      NA
## absences       NA  0.006282370  0.13691739      NA -0.1095387461      NA
## G3             NA  0.103455647 -0.15955029      NA -0.0814071098      NA
##                   Mjob        reason   traveltime   studytime failures
## school              NA            NA           NA          NA       NA
## sex        -0.10277391  0.0098218168  0.020908545  0.18717009       NA
## age         0.08962071 -0.0134006027  0.115165279 -0.05072535       NA
## address             NA            NA           NA          NA       NA
## famsize     0.09157223  0.0009812261 -0.057708148  0.03412913       NA
## Pstatus             NA            NA           NA          NA       NA
## Mjob        1.00000000 -0.0486185847  0.105091345  0.03289251       NA
## reason     -0.04861858  1.0000000000  0.083587795  0.10945933       NA
## traveltime  0.10509135  0.0835877947  1.000000000  0.01494850       NA
## studytime   0.03289251  0.1094593263  0.014948499  1.00000000       NA
## failures            NA            NA           NA          NA        1
## schoolsup           NA            NA           NA          NA       NA
## famsup      0.14047224  0.0762296929  0.006122867  0.08144172       NA
## paid        0.16676148  0.0867430657  0.040906720  0.12737933       NA
## activities  0.12072516  0.0207472217  0.007557387  0.06199645       NA
## nursery             NA            NA           NA          NA       NA
## higher              NA            NA           NA          NA       NA
## internet            NA            NA           NA          NA       NA
## romantic   -0.04106315 -0.0050504430  0.029087153  0.13600710       NA
## famrel      0.02556143  0.0572004768 -0.024250888 -0.04504248       NA
## freetime   -0.05205738  0.0352369814 -0.046113188  0.07137651       NA
## goout      -0.01377800  0.0170752616 -0.037824830  0.01081643       NA
## Dalc       -0.05968205  0.0297294798  0.038855019  0.11420202       NA
## Walc       -0.01954561  0.0602322517  0.030763971  0.12417235       NA
## health     -0.04147519  0.0700755406 -0.019762406  0.05172415       NA
## absences   -0.04718168 -0.0302451714 -0.052442404 -0.02821618       NA
## G3         -0.14572680 -0.0085015237 -0.105385244  0.03982203       NA
##            schoolsup       famsup         paid    activities nursery higher
## school            NA           NA           NA            NA      NA     NA
## sex               NA  0.151622617  0.129125619 -0.0998331770      NA     NA
## age               NA  0.130218615  0.026814641  0.0944094170      NA     NA
## address           NA           NA           NA            NA      NA     NA
## famsize           NA -0.112893067 -0.013882162 -0.0001131763      NA     NA
## Pstatus           NA           NA           NA            NA      NA     NA
## Mjob              NA  0.140472242  0.166761480  0.1207251610      NA     NA
## reason            NA  0.076229693  0.086743066  0.0207472217      NA     NA
## traveltime        NA  0.006122867  0.040906720  0.0075573873      NA     NA
## studytime         NA  0.081441723  0.127379327  0.0619964523      NA     NA
## failures          NA           NA           NA            NA      NA     NA
## schoolsup          1           NA           NA            NA      NA     NA
## famsup            NA  1.000000000  0.293184339 -0.0015001081      NA     NA
## paid              NA  0.293184339  1.000000000 -0.0213823757      NA     NA
## activities        NA -0.001500108 -0.021382376  1.0000000000      NA     NA
## nursery           NA           NA           NA            NA       1     NA
## higher            NA           NA           NA            NA      NA      1
## internet          NA           NA           NA            NA      NA     NA
## romantic          NA  0.012439892  0.005535859  0.0196505434      NA     NA
## famrel            NA  0.025926895  0.036491388 -0.0350101150      NA     NA
## freetime          NA  0.025196125  0.070312152 -0.0532456586      NA     NA
## goout             NA  0.015631443 -0.010493268 -0.0460876858      NA     NA
## Dalc              NA  0.083896122  0.013793327  0.0777486766      NA     NA
## Walc              NA  0.086687935 -0.060453643  0.0374766962      NA     NA
## health            NA -0.029297452  0.078132403 -0.0239228148      NA     NA
## absences          NA -0.024280907 -0.058049073 -0.0526057549      NA     NA
## G3                NA  0.039157145 -0.101996241 -0.0160997013      NA     NA
##            internet     romantic       famrel     freetime        goout
## school           NA           NA           NA           NA           NA
## sex              NA  0.102023009  0.069375471  0.207351483  0.075897396
## age              NA -0.153163415  0.032370369  0.026636252  0.111475815
## address          NA           NA           NA           NA           NA
## famsize          NA  0.034394816  0.006922336 -0.013121698 -0.023064384
## Pstatus          NA           NA           NA           NA           NA
## Mjob             NA -0.041063155  0.025561429 -0.052057378 -0.013777997
## reason           NA -0.005050443  0.057200477  0.035236981  0.017075262
## traveltime       NA  0.029087153 -0.024250888 -0.046113188 -0.037824830
## studytime        NA  0.136007100 -0.045042484  0.071376514  0.010816431
## failures         NA           NA           NA           NA           NA
## schoolsup        NA           NA           NA           NA           NA
## famsup           NA  0.012439892  0.025926895  0.025196125  0.015631443
## paid             NA  0.005535859  0.036491388  0.070312152 -0.010493268
## activities       NA  0.019650543 -0.035010115 -0.053245659 -0.046087686
## nursery          NA           NA           NA           NA           NA
## higher           NA           NA           NA           NA           NA
## internet          1           NA           NA           NA           NA
## romantic         NA  1.000000000  0.045505206  0.032679335 -0.007869926
## famrel           NA  0.045505206  1.000000000  0.093644457  0.087047617
## freetime         NA  0.032679335  0.093644457  1.000000000  0.236439409
## goout            NA -0.007869926  0.087047617  0.236439409  1.000000000
## Dalc             NA  0.067951491 -0.065351965  0.130440135  0.210311226
## Walc             NA  0.010140952 -0.093566930  0.125357672  0.420385745
## health           NA -0.026342317  0.043777985  0.065140450 -0.009577254
## absences         NA -0.059041770 -0.023375472  0.015057208  0.131320084
## G3               NA  0.129969950  0.057689643 -0.003327928 -0.132791474
##                   Dalc        Walc       health    absences           G3
## school              NA          NA           NA          NA           NA
## sex         0.22499694  0.27419377  0.143588173  0.00628237  0.103455647
## age         0.09032711  0.09719130 -0.043727578  0.13691739 -0.159550295
## address             NA          NA           NA          NA           NA
## famsize    -0.14859457 -0.10342501  0.028991659 -0.10953875 -0.081407110
## Pstatus             NA          NA           NA          NA           NA
## Mjob       -0.05968205 -0.01954561 -0.041475191 -0.04718168 -0.145726800
## reason      0.02972948  0.06023225  0.070075541 -0.03024517 -0.008501524
## traveltime  0.03885502  0.03076397 -0.019762406 -0.05244240 -0.105385244
## studytime   0.11420202  0.12417235  0.051724153 -0.02821618  0.039822028
## failures            NA          NA           NA          NA           NA
## schoolsup           NA          NA           NA          NA           NA
## famsup      0.08389612  0.08668793 -0.029297452 -0.02428091  0.039157145
## paid        0.01379333 -0.06045364  0.078132403 -0.05804907 -0.101996241
## activities  0.07774868  0.03747670 -0.023922815 -0.05260575 -0.016099701
## nursery             NA          NA           NA          NA           NA
## higher              NA          NA           NA          NA           NA
## internet            NA          NA           NA          NA           NA
## romantic    0.06795149  0.01014095 -0.026342317 -0.05904177  0.129969950
## famrel     -0.06535196 -0.09356693  0.043777985 -0.02337547  0.057689643
## freetime    0.13044013  0.12535767  0.065140450  0.01505721 -0.003327928
## goout       0.21031123  0.42038575 -0.009577254  0.13132008 -0.132791474
## Dalc        1.00000000  0.53544088  0.074332685  0.05337669 -0.078730632
## Walc        0.53544088  1.00000000  0.092476317  0.17807684 -0.051939324
## health      0.07433269  0.09247632  1.000000000 -0.04993949 -0.061334605
## absences    0.05337669  0.17807684 -0.049939495  1.00000000  0.121485428
## G3         -0.07873063 -0.05193932 -0.061334605  0.12148543  1.000000000

Vemos que hay algunas filas en la matriz de correlación con el valor ‘NA’ correspondientes a las variables school, address, Pstatus, failures, schoolsup, nursery, higher e internet. Esto se debe a que todos los registros de estas variables tienen el mismo valor (en este caso esto se debe al tratamiento de los outliers), haciendo que su desviación típica sea igual a 0 y, por lo tanto, no podemos calcular su correlación con otras variables (ya que sería dividir por 0, lo cual no es posible). Por tanto, debido a que disponemos de muchas más variables, se ha decidido eliminar estas variables.

# The variables 'school', 'address', 'Pstatus', 'failures', 'schoolsup', 'nursery', 'higher' and 'internet' are eliminated
data_xlsx_reduced<-data_xlsx[, -c(1, 4, 6, 11, 12, 16, 17, 18)]

# The first three records in the database are displayed
head(data_xlsx_reduced, n=10)

En cualquier caso, para hacer el análisis de la correlación hemos trabajado con la matriz de correlaciones de todas las variables. Aquí hay un problema porque estamos mezclando variables por un lado continuas o discretas, con variables categóricas que han sido codificadas (sexo, school, address, etc.). Es un problema porque por defecto estamos obteniendo las correlaciones de Pearson que trata todas las variables como cuantitativas y no lo son, llevando así a engaño. Para variables categóricas habría que cambiar el parámetro de la función de R para que obtuviera la correlación de Spearman. Pero claro, ya no podríamos hacer una reducción de la dimensión con todas a la vez. Por un lado lo haríamos con las cuantitativas que es para lo que funciona el ACP por su construcción y por otro lado investigaríamos como hacerlo con variables categóricas (lo cual no hemos visto). Este comentario es solo para tenerlo en cuenta de cara al futuro. Para esta práctica, el procedimiento seguido mezclando todas las variables es suficiente. En conclusión, recordar para trabajos futuros que la reducción de la dimensión la hacéis con las variables del mismo tipo (cuantitativas).

De acuerdo con los resultados numéricos a continuación, se observa que los datos no están correlacionados ni a nivel muestral (ver matriz de correlación) ni a nivel poblacional (la prueba de esfericidad de Bartlett no es significativa).

###############################
# Correlation at sample level #
###############################

# Are the variables correlated at sample level?
correlation_matrix<-cor(data_xlsx_reduced)
correlation_matrix
##                     sex         age       famsize        Mjob        reason
## sex         1.000000000 -0.04064933 -0.0898618018 -0.10277391  0.0098218168
## age        -0.040649328  1.00000000 -0.0455885843  0.08962071 -0.0134006027
## famsize    -0.089861802 -0.04558858  1.0000000000  0.09157223  0.0009812261
## Mjob       -0.102773909  0.08962071  0.0915722328  1.00000000 -0.0486185847
## reason      0.009821817 -0.01340060  0.0009812261 -0.04861858  1.0000000000
## traveltime  0.020908545  0.11516528 -0.0577081482  0.10509135  0.0835877947
## studytime   0.187170093 -0.05072535  0.0341291336  0.03289251  0.1094593263
## famsup      0.151622617  0.13021861 -0.1128930671  0.14047224  0.0762296929
## paid        0.129125619  0.02681464 -0.0138821617  0.16676148  0.0867430657
## activities -0.099833177  0.09440942 -0.0001131763  0.12072516  0.0207472217
## romantic    0.102023009 -0.15316341  0.0343948162 -0.04106315 -0.0050504430
## famrel      0.069375471  0.03237037  0.0069223355  0.02556143  0.0572004768
## freetime    0.207351483  0.02663625 -0.0131216979 -0.05205738  0.0352369814
## goout       0.075897396  0.11147582 -0.0230643843 -0.01377800  0.0170752616
## Dalc        0.224996943  0.09032711 -0.1485945736 -0.05968205  0.0297294798
## Walc        0.274193767  0.09719130 -0.1034250076 -0.01954561  0.0602322517
## health      0.143588173 -0.04372758  0.0289916591 -0.04147519  0.0700755406
## absences    0.006282370  0.13691739 -0.1095387461 -0.04718168 -0.0302451714
## G3          0.103455647 -0.15955029 -0.0814071098 -0.14572680 -0.0085015237
##              traveltime   studytime       famsup         paid    activities
## sex         0.020908545  0.18717009  0.151622617  0.129125619 -0.0998331770
## age         0.115165279 -0.05072535  0.130218615  0.026814641  0.0944094170
## famsize    -0.057708148  0.03412913 -0.112893067 -0.013882162 -0.0001131763
## Mjob        0.105091345  0.03289251  0.140472242  0.166761480  0.1207251610
## reason      0.083587795  0.10945933  0.076229693  0.086743066  0.0207472217
## traveltime  1.000000000  0.01494850  0.006122867  0.040906720  0.0075573873
## studytime   0.014948499  1.00000000  0.081441723  0.127379327  0.0619964523
## famsup      0.006122867  0.08144172  1.000000000  0.293184339 -0.0015001081
## paid        0.040906720  0.12737933  0.293184339  1.000000000 -0.0213823757
## activities  0.007557387  0.06199645 -0.001500108 -0.021382376  1.0000000000
## romantic    0.029087153  0.13600710  0.012439892  0.005535859  0.0196505434
## famrel     -0.024250888 -0.04504248  0.025926895  0.036491388 -0.0350101150
## freetime   -0.046113188  0.07137651  0.025196125  0.070312152 -0.0532456586
## goout      -0.037824830  0.01081643  0.015631443 -0.010493268 -0.0460876858
## Dalc        0.038855019  0.11420202  0.083896122  0.013793327  0.0777486766
## Walc        0.030763971  0.12417235  0.086687935 -0.060453643  0.0374766962
## health     -0.019762406  0.05172415 -0.029297452  0.078132403 -0.0239228148
## absences   -0.052442404 -0.02821618 -0.024280907 -0.058049073 -0.0526057549
## G3         -0.105385244  0.03982203  0.039157145 -0.101996241 -0.0160997013
##                romantic       famrel     freetime        goout        Dalc
## sex         0.102023009  0.069375471  0.207351483  0.075897396  0.22499694
## age        -0.153163415  0.032370369  0.026636252  0.111475815  0.09032711
## famsize     0.034394816  0.006922336 -0.013121698 -0.023064384 -0.14859457
## Mjob       -0.041063155  0.025561429 -0.052057378 -0.013777997 -0.05968205
## reason     -0.005050443  0.057200477  0.035236981  0.017075262  0.02972948
## traveltime  0.029087153 -0.024250888 -0.046113188 -0.037824830  0.03885502
## studytime   0.136007100 -0.045042484  0.071376514  0.010816431  0.11420202
## famsup      0.012439892  0.025926895  0.025196125  0.015631443  0.08389612
## paid        0.005535859  0.036491388  0.070312152 -0.010493268  0.01379333
## activities  0.019650543 -0.035010115 -0.053245659 -0.046087686  0.07774868
## romantic    1.000000000  0.045505206  0.032679335 -0.007869926  0.06795149
## famrel      0.045505206  1.000000000  0.093644457  0.087047617 -0.06535196
## freetime    0.032679335  0.093644457  1.000000000  0.236439409  0.13044013
## goout      -0.007869926  0.087047617  0.236439409  1.000000000  0.21031123
## Dalc        0.067951491 -0.065351965  0.130440135  0.210311226  1.00000000
## Walc        0.010140952 -0.093566930  0.125357672  0.420385745  0.53544088
## health     -0.026342317  0.043777985  0.065140450 -0.009577254  0.07433269
## absences   -0.059041770 -0.023375472  0.015057208  0.131320084  0.05337669
## G3          0.129969950  0.057689643 -0.003327928 -0.132791474 -0.07873063
##                   Walc       health    absences           G3
## sex         0.27419377  0.143588173  0.00628237  0.103455647
## age         0.09719130 -0.043727578  0.13691739 -0.159550295
## famsize    -0.10342501  0.028991659 -0.10953875 -0.081407110
## Mjob       -0.01954561 -0.041475191 -0.04718168 -0.145726800
## reason      0.06023225  0.070075541 -0.03024517 -0.008501524
## traveltime  0.03076397 -0.019762406 -0.05244240 -0.105385244
## studytime   0.12417235  0.051724153 -0.02821618  0.039822028
## famsup      0.08668793 -0.029297452 -0.02428091  0.039157145
## paid       -0.06045364  0.078132403 -0.05804907 -0.101996241
## activities  0.03747670 -0.023922815 -0.05260575 -0.016099701
## romantic    0.01014095 -0.026342317 -0.05904177  0.129969950
## famrel     -0.09356693  0.043777985 -0.02337547  0.057689643
## freetime    0.12535767  0.065140450  0.01505721 -0.003327928
## goout       0.42038575 -0.009577254  0.13132008 -0.132791474
## Dalc        0.53544088  0.074332685  0.05337669 -0.078730632
## Walc        1.00000000  0.092476317  0.17807684 -0.051939324
## health      0.09247632  1.000000000 -0.04993949 -0.061334605
## absences    0.17807684 -0.049939495  1.00000000  0.121485428
## G3         -0.05193932 -0.061334605  0.12148543  1.000000000
det(correlation_matrix)
## [1] 0.1963095
###################################
# Correlation at population level #
###################################

# Bartlett's sphericity test:
# This test checks whether the correlations are significantly different from 0
# The null hypothesis is H_0; det(R)=1 means the variables are uncorrelated 
# R denotes the correlation matrix
# cortest.bartlett function in the package pysch performs this test
# This function works with standardized data.

# Standardization
data_xlsx_scale<-scale(data_xlsx_reduced)

# Bartlett's sphericity test
cortest.bartlett(cor(data_xlsx_scale))
## $chisq
## [1] 149.5104
## 
## $p.value
## [1] 0.880634
## 
## $df
## [1] 171

A la vista de los resultados del test (p-value > 0.001), no se puede rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, no tendría mucho sentido realmente plantearse una reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF). Sin embargo, vamos a realizarlo con el fin de mostrar cómo se haría.

3.2 Realización del ACP

El siguiente código realiza el ACP, obteniendo los vectores propios que generan cada componente, así como sus valores propios, que corresponden a la varianza de cada una.

# The 'prcomp' function in the base R package performs this analysis
# Parameters 'scale' and 'center' are set to TRUE to consider standardized data
PCA<-prcomp(data_xlsx_reduced, scale=T, center=T)

# The field 'rotation' of the 'PCA' object is a matrix 
# Its columns are the coefficients of the principal components
# Indicates the weight of each variable in the corresponding principal component
PCA$rotation
##                     PC1          PC2          PC3          PC4          PC5
## sex         0.357988987 -0.138411478  0.322697872 -0.054267662  0.067068174
## age         0.131267998  0.339957291 -0.331365990 -0.176090198  0.145042062
## famsize    -0.161590831  0.071228092  0.069841422  0.498471643  0.008965111
## Mjob       -0.044001488  0.492619024  0.007931548 -0.002584637 -0.034385460
## reason      0.094556511  0.101828074  0.186479728  0.076449751 -0.014264218
## traveltime  0.038760178  0.271918073 -0.015201624 -0.029883779 -0.201107849
## studytime   0.194298544  0.043651596  0.346186243  0.055168073 -0.299681437
## famsup      0.182311516  0.318048173  0.250930379 -0.422734663  0.118971288
## paid        0.104874161  0.392601699  0.366688720 -0.111221050  0.194371680
## activities -0.001007259  0.198060733 -0.082199948 -0.009605777 -0.425364174
## romantic    0.054747386 -0.182091711  0.321552812  0.017569240 -0.268945880
## famrel      0.008972162 -0.007201713  0.152799749  0.037058002  0.483625472
## freetime    0.270209598 -0.087769965  0.102905512  0.197956338  0.343553056
## goout       0.366548664 -0.025734817 -0.255907511  0.205602698  0.234858149
## Dalc        0.465023159 -0.022189337 -0.104584992  0.036943252 -0.268693347
## Walc        0.519812584 -0.059801406 -0.197971823  0.062367826 -0.178380077
## health      0.123665610 -0.012696051  0.178673612  0.349364636  0.083594720
## absences    0.139648827 -0.166723530 -0.298530527 -0.351770959  0.135539274
## G3         -0.049366010 -0.403967956  0.219294237 -0.427111617 -0.032997838
##                    PC6         PC7         PC8          PC9        PC10
## sex        -0.10391169 -0.10839454  0.16365183 -0.139257379  0.20557884
## age         0.02303033  0.15411088 -0.07849430 -0.088680141  0.22237639
## famsize     0.31314158 -0.10650027 -0.10643092  0.301916325  0.25727557
## Mjob        0.33116460 -0.11367398  0.06714274 -0.011128482  0.24163018
## reason     -0.28271349  0.52853728 -0.45083302  0.395241591 -0.25834592
## traveltime -0.28597132  0.48146496  0.49939418  0.102358812  0.29838847
## studytime   0.14566345 -0.03532613 -0.19240446  0.298003384  0.20817459
## famsup      0.04890066 -0.15759837 -0.04364686  0.004294231 -0.21485205
## paid       -0.04884339 -0.22592945 -0.01539739  0.099576591 -0.05774483
## activities  0.29229668  0.16916938 -0.43267571 -0.456748314 -0.06963954
## romantic    0.37116072  0.22067039  0.32392297 -0.002982180  0.04013603
## famrel      0.23783650  0.48463424 -0.04143655 -0.345570359  0.13434361
## freetime    0.18335791  0.05011800  0.04997650 -0.026074757 -0.15056639
## goout       0.26727354  0.06983414  0.05217576  0.171907752 -0.17261233
## Dalc       -0.05127171 -0.03837628  0.05303724 -0.175134786 -0.11882933
## Walc        0.01527602 -0.04995084 -0.02075521  0.044733675  0.02095125
## health     -0.41945363 -0.13236526 -0.27465340 -0.343906254  0.42252253
## absences    0.07709932 -0.02275517 -0.22455581  0.325558847  0.48661653
## G3          0.15387047  0.10182366 -0.17674737 -0.028642312  0.16806270
##                    PC11         PC12         PC13         PC14        PC15
## sex        -0.143978916  0.197058675 -0.263363032 -0.002061750 -0.30748633
## age        -0.324515844  0.487108748  0.246267697 -0.055017486  0.22599694
## famsize     0.035308463  0.488004212 -0.100712430  0.301302539 -0.24592511
## Mjob        0.242299180 -0.380084492 -0.340510269  0.140047267  0.23154673
## reason      0.136849347 -0.006144228 -0.028546872  0.161250872  0.01531048
## traveltime -0.128723700 -0.096826268 -0.131715638  0.175449808 -0.08922372
## studytime  -0.306801210 -0.053734447 -0.089034545 -0.600946157  0.23786877
## famsup      0.169697930  0.312720452 -0.037366643  0.189046850  0.28369633
## paid        0.013212497 -0.190954934  0.330258017 -0.002415829 -0.48699543
## activities -0.195063366 -0.123660716  0.051135496  0.131056496 -0.27602362
## romantic    0.205147170  0.046890289  0.599315081  0.119626838  0.17559823
## famrel      0.253285005  0.053805328 -0.121153360 -0.360065071 -0.15430665
## freetime   -0.608901579 -0.273425928  0.012434492  0.362500709  0.14683036
## goout       0.173901756 -0.130000373  0.008027874 -0.067047736  0.10288304
## Dalc        0.133241456  0.105844736  0.006299095  0.002167397 -0.16751420
## Walc        0.240014167  0.043755572 -0.175019668  0.064239619 -0.01442977
## health      0.177947027 -0.109242105  0.214537433  0.158185734  0.33059549
## absences    0.057199275 -0.238459491  0.267427336  0.061761636 -0.21955060
## G3         -0.004378558  0.027574136 -0.287120761  0.320784978  0.09547792
##                   PC16        PC17        PC18        PC19
## sex         0.04393326 -0.59012512  0.15568682  0.17953965
## age         0.08134033  0.01016004  0.39072602 -0.02951198
## famsize    -0.03244157  0.13339112 -0.12771368  0.03680353
## Mjob        0.32382919 -0.11862328  0.22129067  0.09097057
## reason      0.21928319 -0.18524533  0.14898620  0.06297112
## traveltime -0.27960227  0.15174886 -0.19130319  0.02558064
## studytime  -0.07042271  0.13531241 -0.09935796  0.01473493
## famsup     -0.15542534 -0.04212068 -0.51168242  0.05009488
## paid       -0.15142989  0.26956896  0.29247492 -0.13899018
## activities -0.24644475 -0.17754971 -0.12756713  0.02002253
## romantic    0.10271605 -0.16220615  0.08151722 -0.05188973
## famrel      0.10093130  0.14617421 -0.17382129 -0.11696431
## freetime    0.20357093  0.09979195 -0.16600066 -0.10908413
## goout      -0.57289816 -0.07682474  0.17385188  0.37020706
## Dalc        0.40773927  0.46952353 -0.04379517  0.44874450
## Walc       -0.03853755  0.01675788  0.04086588 -0.74013052
## health     -0.16081574  0.05391938 -0.03258438  0.06412728
## absences    0.13395249 -0.11290448 -0.31340469  0.09598673
## G3         -0.21108835  0.37180623  0.35777495  0.04881103
# Standard deviations of each principal component
PCA$sdev
##  [1] 1.4966008 1.2622886 1.2484457 1.1234160 1.1212370 1.0524091 1.0247914
##  [8] 1.0062424 0.9838951 0.9547318 0.9281544 0.8978759 0.8823035 0.8702234
## [15] 0.8226476 0.8007155 0.7737687 0.7478517 0.5998485

Cada componente principal se obtiene de forma sencilla como una combinación lineal de todas las variables con los coeficientes indicados por las columnas de la matriz de rotación.

3.3 Tasa de variación explicada

# The function 'summary' applied to the 'PCA' object provides relevant information
# - Standard deviations of each principal component
# - Proportion of variance explained and cummulative variance
summary(PCA)
## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.4966 1.26229 1.24845 1.12342 1.12124 1.05241 1.02479
## Proportion of Variance 0.1179 0.08386 0.08203 0.06642 0.06617 0.05829 0.05527
## Cumulative Proportion  0.1179 0.20175 0.28378 0.35020 0.41637 0.47466 0.52994
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     1.00624 0.98390 0.95473 0.92815 0.89788 0.88230 0.87022
## Proportion of Variance 0.05329 0.05095 0.04797 0.04534 0.04243 0.04097 0.03986
## Cumulative Proportion  0.58323 0.63418 0.68215 0.72749 0.76992 0.81089 0.85075
##                           PC15    PC16    PC17    PC18    PC19
## Standard deviation     0.82265 0.80072 0.77377 0.74785 0.59985
## Proportion of Variance 0.03562 0.03374 0.03151 0.02944 0.01894
## Cumulative Proportion  0.88637 0.92011 0.95163 0.98106 1.00000

Los siguientes gráficos ilustran el comportamiento de la varianza explicada por cada componente principal, así como el comportamiento de la varianza explicada acumulada.

# The following graph shows the proportion of explained variance
Explained_variance <- PCA$sdev^2 / sum(PCA$sdev^2)

p13<-ggplot(data = data.frame(Explained_variance, pc = 1:ncol(data_xlsx_reduced)),
  aes(x = pc, y = Explained_variance, fill=Explained_variance )) +
  geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
  labs(x = "Principal component", y= "Proportion of variance")

# The following graph shows the proportion of cumulative explained variance
Cummulative_variance<-cumsum(Explained_variance)

p14<-ggplot( data = data.frame(Cummulative_variance, pc = 1:ncol(data_xlsx_reduced)),
  aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) +
  geom_col(width = 0.5) +  scale_y_continuous(limits = c(0,1)) + theme_bw() +
  labs(x = "Principal component",
       y = "Proportion of cumulative variance")

p13

p14

3.4 Número de componentes principales apropiado

Para elegir el número adecuado de componentes principales, existen diferentes métodos. Nosotros vamos a utilizar la Regla de Abdi et al. (2010): Se promedian las varianzas explicadas por las componentes principales y se seleccionan aquellas cuya proporción de varianza explicada supera la media.

#######################
# Rule of Abdi et al. #
#######################

# Variances
PCA$sdev^2
##  [1] 2.2398140 1.5933726 1.5586166 1.2620636 1.2571724 1.1075649 1.0501974
##  [8] 1.0125238 0.9680496 0.9115128 0.8614706 0.8061811 0.7784594 0.7572887
## [15] 0.6767490 0.6411453 0.5987180 0.5592821 0.3598182
# Average of variances
mean(PCA$sdev^2)
## [1] 1

A la vista de los resultados, en este caso, se eligen ocho direcciones principales que, tal y como se puede ver, acumulan alrededor del 60% de varianza explicada.

3.5 Salidas gráficas de interés

Los siguientes gráficos muestran la comparativa entre distintas componentes principales mediante una proyección al plano sobre cada dos componentes. En esta representación se aprecia cuáles de las variables originales tienen mayor o menor peso en cada una de las componentes enfrentadas.

# These graphical outputs show the projection of the variables in two dimensions
# Display the weight of the variable in the direction of the principal component 
p15<-fviz_pca_var(PCA, repel=TRUE, col.var="cos2",
                 legend.title="Distance", title="Variables")+theme_bw()

p16<-fviz_pca_var(PCA, axes=c(1,3), repel=TRUE, col.var="cos2",
                 legend.title="Distance", title="Variables")+theme_bw()

p17<-fviz_pca_var(PCA, axes=c(2,3), repel=TRUE, col.var="cos2",
                 legend.title="Distance", title="Variables")+theme_bw()

# Displaying graphics
p15

p16

p17

Es posible también representar las observaciones de los objetos junto con las componentes principales mediante la orden contrib() de la función fviz_pca_ind() anterior, así como identificar con colores aquellas observaciones que mayor varianza explican.

# It is also possible to represent the observations
# As well as identify with colors those observations that explain the greatest 
# variance of the principal components
p18<-fviz_pca_ind(PCA,col.ind = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

p19<-fviz_pca_ind(PCA,axes=c(1,3),col.ind = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

p20<-fviz_pca_ind(PCA,axes=c(2,3),col.ind = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

# Displaying graphics
p18

p19

p20

Finalmente, se puede obtener una representación conjunta de variables y observaciones que relaciona visualmente las posibles relaciones entre las observaciones, las contribuciones de los individuos a las varianzas y el peso de las variables en cada componente principal.

# Joint representation of variables and observations
# Relates the possible relationships between the contributions of the records
# to the variances of the components and the weight of the variables in each 
# principal component

p21<-fviz_pca(PCA,alpha.ind ="contrib", col.var = "cos2",
         col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE, legend.title="Distancia")+theme_bw()

p22<-fviz_pca(PCA,axes=c(1,3),alpha.ind ="contrib", 
         col.var = "cos2",col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE, legend.title="Distancia")+theme_bw()

p23<-fviz_pca(PCA,axes=c(2,3),alpha.ind ="contrib", 
         col.var = "cos2",col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE, legend.title="Distancia")+theme_bw()

# Displaying graphics
p21

p22

p23

Finalmente, dado que el objeto de este estudio es reducir la dimensión del conjunto de datos, es posible obtener las coordenadas de los datos originales en el nuevo sistema de referencia. De hecho, se almacenan desde que usamos la función prcomp() para crear la variable PCA.

head(PCA$x)
##             PC1        PC2         PC3        PC4        PC5        PC6
## [1,] -0.5428656  2.2596133 -0.42483568 -0.3808055  0.1973962  0.5663425
## [2,] -1.4738393  0.9734234 -0.03380767  0.7280407  0.6325970  1.0126359
## [3,]  0.3584064 -0.4617297 -0.43632964 -1.6441159 -1.2831378 -0.2222366
## [4,] -2.7481266 -1.8532500 -0.57242625  0.3913111 -0.6674298 -1.6632617
## [5,] -1.8208294 -0.3271368 -0.59114102  0.8691652 -0.8624869  0.8546701
## [6,] -0.0997650 -2.5648834  0.38175788 -0.9531831  1.4412018 -0.5339585
##             PC7         PC8          PC9       PC10      PC11       PC12
## [1,]  0.9718671 -0.04616334  0.597973268 -0.2713659 0.2946118  0.3563613
## [2,]  1.0458847 -0.73175801 -0.358337362 -0.4003332 0.6240971 -0.1919391
## [3,]  0.8209019 -1.31891932 -0.048872389 -1.2455403 1.6704745 -1.0877667
## [4,] -2.0681437 -0.21915098 -0.006630797  0.2687390 0.3567702  0.3884771
## [5,] -0.7398849 -0.20685083 -1.425692474  0.9659082 0.6450230 -0.5232906
## [6,]  0.5215257  0.11646362 -1.123298254  1.2918069 0.2256947 -0.9883929
##            PC13        PC14       PC15       PC16       PC17        PC18
## [1,]  1.4946749  0.89213985 -0.2480674 -0.8491035 -0.4345776 -0.45757939
## [2,]  1.2684351 -0.26645691 -1.0759564  0.4760184 -0.2633952  0.33962636
## [3,]  0.4009926  0.80278031  0.2609793  1.6434488 -0.5685151 -1.55488057
## [4,] -0.6810208  0.17717504  0.4006045 -0.7550371  0.8415472 -0.06199315
## [5,]  0.3234729  0.71503870  0.5434042  0.4081917 -0.2519953 -0.30382251
## [6,]  0.2993773  0.06754642  0.4105511  1.1864022 -0.8217480 -0.16359497
##             PC19
## [1,]  0.65503362
## [2,] -0.02653478
## [3,]  0.03444780
## [4,]  0.40193837
## [5,] -0.38369841
## [6,] -0.31557997

Recordemos que no hemos podido rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente.

4 Análisis Factorial (AF)

4.1 Correlación

Es necesario comprobar que las variables no son independientes. Ya vimos en el apartado 3, al comprobar la correlación de las variables para el Análisis de Componentes Principales, que las variables sí son independientes tras calcular la matriz de correlación y realizar el test de Bartlett. Esto nos dice, volvemos a repetirlo, qque la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente. Sin embargo, igual que hicimos en el Análisis de Componentes Principales, vamos a realizarlo con el fin de mostrar cómo se haría.

Caben destacar también los siguientes resultados gráficos que proporcionan una idea intuitiva de la correlación entre las variables:

# Polychoric correlation matrix
# Remember that package 'polycor' is required for 'hetcor' function
poly_cor<-hetcor(data_xlsx_reduced)$correlations

# Remember that package 'ggcoorplot' is required for 'ggcorrplot' function
ggcorrplot(poly_cor, type="lower", hc.order=T)

# Another interesting visual representation is the following
# Remember that package 'corrplot' is required for 'corrplot' function
corrplot(cor(data_xlsx_reduced), order = "hclust", tl.col='black', tl.cex=1)

# Or this other option is also very visual
# Remember that package 'corrr' is required for 'rplot' function
data_xlsx_reduced_correlations <- correlate(data_xlsx_reduced)
rplot(data_xlsx_reduced_correlations, legend = TRUE, colours = c("firebrick1", "black","darkcyan"), print_cor = TRUE) 

Tal y como era de esperar, podemos observar que las variables apenas están correladas entre sí (lo cual corrobora los resultados numéricos de la matriz de correlación y el test de Bartlett que obtuvimos anteriormente).

4.2 Diferentes métodos de estimación

Para extraer los factores hay que escoger un método de estimación. Para ello, se utiliza la función fa(), que permite implementar hasta 6 métodos distintos.

A continuación, se comparan las salidas com el método del factor principal y con el método de máxima verosimilitud. Se calcula 1 factor, ya que es lo que se ve intuitivamente en los mapas de correlaciones anteriores.

### Test of two models with three factors
modelo1<-fa(poly_cor,
            nfactors = 1,
            rotate = "none",
            fm="mle") # likelihood method

modelo2<-fa(poly_cor,
            nfactors = 1,
            rotate = "none",
            fm="minres") # minimal residual model

# Outputs of these models: factorial matrices, etc.

print("Modelo 1: mle")
## [1] "Modelo 1: mle"
modelo1$loadings
## 
## Loadings:
##            ML1   
## sex         0.320
## age         0.123
## famsize    -0.136
## Mjob             
## reason           
## traveltime       
## studytime   0.151
## famsup      0.112
## paid             
## activities       
## romantic         
## famrel           
## freetime    0.185
## goout       0.455
## Dalc        0.603
## Walc        0.884
## health      0.105
## absences    0.185
## G3               
## 
##                  ML1
## SS loadings    1.623
## Proportion Var 0.085
print("Modelo 2: minres")
## [1] "Modelo 2: minres"
modelo2$loadings
## 
## Loadings:
##            MR1   
## sex         0.376
## age         0.134
## famsize    -0.163
## Mjob             
## reason           
## traveltime       
## studytime   0.185
## famsup      0.161
## paid             
## activities       
## romantic         
## famrel           
## freetime    0.264
## goout       0.432
## Dalc        0.615
## Walc        0.809
## health      0.119
## absences    0.156
## G3               
## 
##                  MR1
## SS loadings    1.596
## Proportion Var 0.084

Se comparan las comunalidades de ambos modelos para ver qué proporción de la varianza es explicada por los factores latentes.

# Comparing communalities
sort(modelo1$communality, decreasing = T)->c1
sort(modelo2$communality, decreasing = T)->c2
head(cbind(c1, c2))
##                  c1         c2
## Walc     0.78204940 0.65515285
## Dalc     0.36417822 0.37848395
## goout    0.20704015 0.18671636
## sex      0.10217693 0.14110847
## freetime 0.03418673 0.06982291
## absences 0.03415778 0.03429621

También se comparan las unicidades, esto es, la proporción de varianza que no ha sido explicada por el factor (1-comunalidad):

# Comparison of uniqueness, that is, the proportion of variance that has not 
# been explained by the factor (1-communality)
sort(modelo1$uniquenesses, decreasing = T)->u1
sort(modelo2$uniquenesses, decreasing = T)->u2
head(cbind(u1,u2), n=10)
##                   u1        u2
## romantic   0.9993603 0.9998923
## paid       0.9993580 0.9997853
## traveltime 0.9988124 0.9986512
## activities 0.9987844 0.9976418
## Mjob       0.9985758 0.9975256
## G3         0.9958093 0.9967005
## reason     0.9954460 0.9956642
## famrel     0.9943776 0.9930515
## health     0.9890089 0.9859061
## famsup     0.9875302 0.9819405

Se observa que con el método de máxima verosimilitud se obtienen valores un poco más altos, luego los factores latentes obtenidos con este método explican mejor la varianza de las variables observadas.

4.3 Número de factores latentes apropiado

Existen diferentes criterios, entre los que destacan el scree plot (Cattel 1966) y el análisis paralelo (Horn 1965).

El método del codo consiste en representar el gráfico de sedimentación, que se obtiene al representar en las ordenadas los valores propios en orden decreciente y en las abscias los números correspondiente a las componentes principales o factores latentes (según se esté haciendo un ACP o un AF). Uniendo los puntos se obtiene una figura que empieza con fuerte pendiente para luego seguir con una ligera inclinación. El método indica que hay que quedarse con el número de factores que haya hasta el “codo” de la gráfica, donde se pasa de la fuerte pendiente a la ligera inclinación.

Se llama gráfico de sedimentación (scree plot en inglés), ya que parece el perfil de una montaña con una gran pendiente hasta llegar a la base, que es una meseta donde se acumulan los guijarros caídos desde la cumbre (donde sedimentan).

Según los siguientes resultados gráficos, 1 se considera el número óptimo de factores latentes (análisis paralelo), coincidiendo con el primer gráfico scree plot.

# Scree plot
scree(poly_cor)

# Parallel analysis
fa.parallel(poly_cor, n.obs=100, fa="fa", fm="ml")

## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

Otra forma de estudiar el número óptimo de factores es con el siguiente test de hipótesis, que contrasta si el número de factores es suficiente o no:

# Remember that package 'stats' is required for 'factanal' function
# This function only performs the mle method
factanal(data_xlsx_reduced, factors=1, rotation="varimax")
## 
## Call:
## factanal(x = data_xlsx_reduced, factors = 1, rotation = "varimax")
## 
## Uniquenesses:
##        sex        age    famsize       Mjob     reason traveltime  studytime 
##      0.898      0.985      0.981      0.999      0.995      0.999      0.977 
##     famsup       paid activities   romantic     famrel   freetime      goout 
##      0.988      0.999      0.999      0.999      0.994      0.966      0.793 
##       Dalc       Walc     health   absences         G3 
##      0.636      0.218      0.989      0.966      0.996 
## 
## Loadings:
##            Factor1
## sex         0.320 
## age         0.123 
## famsize    -0.136 
## Mjob              
## reason            
## traveltime        
## studytime   0.151 
## famsup      0.112 
## paid              
## activities        
## romantic          
## famrel            
## freetime    0.185 
## goout       0.455 
## Dalc        0.603 
## Walc        0.884 
## health      0.105 
## absences    0.185 
## G3                
## 
##                Factor1
## SS loadings      1.623
## Proportion Var   0.085
## 
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 341.42 on 152 degrees of freedom.
## The p-value is 1.36e-16

Así, se acepta la hipótesis de que 1 factor es suficiente.

4.4 Realización del AF

Finalmente, se estima el modelo factorial con 1 factor implementando una rotación de tipo varimax para buscar una interpretación más simple de la realidad.

modelo_varimax<-fa(poly_cor, nfactors=1, rotate="varimax", fa="mle")

# The rotated factorial matrix is shown
print(modelo_varimax$loadings, cut=0) 
## 
## Loadings:
##            MR1   
## sex         0.376
## age         0.134
## famsize    -0.163
## Mjob       -0.050
## reason      0.083
## traveltime  0.037
## studytime   0.185
## famsup      0.161
## paid        0.066
## activities  0.010
## romantic    0.049
## famrel     -0.015
## freetime    0.264
## goout       0.432
## Dalc        0.615
## Walc        0.809
## health      0.119
## absences    0.156
## G3         -0.057
## 
##                  MR1
## SS loadings    1.596
## Proportion Var 0.084

Visualmente se podría hacer el esfuerzo de ver con qué variables correlaciona cada uno de los factores en la matriz factorial, pero esto es muy tedioso, por lo que se utiliza la siguiente representación:

fa.diagram(modelo_varimax)

En el diagrama se observa que el único factor latente que tenemos correlaciona con las variables Walc, Dalc, goout y sex, siendo un resultado muy lógico que era de esperar.

Vemos que solo tenemos un factor latente y que solo correlaciona cuatro variables. Esto se debe a que (lo volvemos a recordar por última vez) no hemos podido rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente.

5 Análisis Discriminante Lineal (ADL)

La base de datos contiene información sobre la nota final del estudiante en la asignatura Matemáticas. Así, sería interesante dar un método de clasificación para la calificación según los demás indicadores.

Para ello, haremos que la variable G3 sea categórica, con dos categorías:

  • Aprobado si la nota dada en la variable G3 es mayor o igual que 10.

  • Suspenso si la nota dada en la variable G3 es menor que 10.

calification = c()

for(k in 1:nrow(data_xlsx)) {
  if(data_xlsx_reduced$G3[k] >= 10){
    calification[k] = "aprobado"
  }
  else
    calification[k] = "suspenso"
}

calification<-as.factor(calification)

# 'G3' now is a cathegorical variable
data_xlsx_G3_cathegorical = data_xlsx_reduced
data_xlsx_G3_cathegorical$G3 = calification

data_xlsx_G3_cathegorical$G3
##   [1] suspenso suspenso aprobado aprobado aprobado aprobado aprobado suspenso
##   [9] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
##  [17] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
##  [25] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
##  [33] aprobado aprobado aprobado suspenso aprobado aprobado aprobado aprobado
##  [41] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
##  [49] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
##  [57] aprobado aprobado suspenso aprobado aprobado aprobado suspenso suspenso
##  [65] aprobado aprobado aprobado suspenso suspenso aprobado aprobado aprobado
##  [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado suspenso
##  [81] aprobado aprobado suspenso aprobado aprobado suspenso suspenso aprobado
##  [89] aprobado suspenso suspenso aprobado suspenso aprobado aprobado aprobado
##  [97] aprobado aprobado aprobado suspenso suspenso aprobado aprobado suspenso
## [105] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
## [113] aprobado aprobado suspenso aprobado aprobado aprobado suspenso aprobado
## [121] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [129] suspenso aprobado suspenso suspenso aprobado aprobado suspenso suspenso
## [137] suspenso suspenso aprobado aprobado suspenso suspenso aprobado aprobado
## [145] suspenso aprobado suspenso aprobado suspenso aprobado suspenso aprobado
## [153] aprobado suspenso aprobado suspenso aprobado aprobado aprobado aprobado
## [161] suspenso suspenso suspenso aprobado suspenso aprobado aprobado aprobado
## [169] suspenso aprobado suspenso aprobado aprobado suspenso suspenso suspenso
## [177] aprobado suspenso suspenso aprobado suspenso aprobado aprobado suspenso
## [185] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [193] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [201] aprobado aprobado aprobado suspenso aprobado suspenso suspenso aprobado
## [209] aprobado suspenso suspenso aprobado aprobado suspenso aprobado aprobado
## [217] suspenso suspenso suspenso aprobado suspenso suspenso aprobado aprobado
## [225] aprobado suspenso aprobado aprobado suspenso aprobado aprobado aprobado
## [233] suspenso aprobado suspenso aprobado aprobado aprobado aprobado suspenso
## [241] aprobado aprobado suspenso aprobado suspenso aprobado aprobado suspenso
## [249] suspenso aprobado suspenso aprobado suspenso suspenso aprobado suspenso
## [257] aprobado aprobado aprobado suspenso aprobado suspenso aprobado suspenso
## [265] suspenso aprobado aprobado aprobado aprobado suspenso suspenso aprobado
## [273] aprobado aprobado aprobado aprobado suspenso suspenso suspenso aprobado
## [281] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [289] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [297] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [305] aprobado aprobado aprobado suspenso aprobado aprobado suspenso aprobado
## [313] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
## [321] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [329] suspenso aprobado suspenso aprobado suspenso suspenso suspenso aprobado
## [337] aprobado suspenso aprobado aprobado aprobado suspenso aprobado suspenso
## [345] aprobado aprobado aprobado suspenso aprobado aprobado suspenso aprobado
## [353] suspenso suspenso aprobado suspenso aprobado aprobado aprobado aprobado
## [361] aprobado aprobado aprobado aprobado aprobado aprobado aprobado suspenso
## [369] aprobado aprobado suspenso aprobado aprobado suspenso aprobado aprobado
## [377] aprobado aprobado aprobado aprobado aprobado suspenso aprobado suspenso
## [385] suspenso aprobado suspenso suspenso suspenso suspenso suspenso aprobado
## [393] suspenso aprobado suspenso
## Levels: aprobado suspenso

Por otro lado, sabemos que el AD se basa en el concepto de distancia de Mahalanobis y se utiliza comúnmente para clasificar casos en grupos. Para realizar un AD, necesitamos al menos una variable categórica que actúe como tu variable dependiente (esta es la variable data_xlsx_G3_cathegorical$G3 en nuestro caso) y varias variables independientes numéricas que se usarán para predecir el grupo de clasificación.

Las variables categóricas binarias también pueden ser incluidas como variables independientes en el AD si se codifican adecuadamente, aunque su interpretación puede ser menos clara que con las variables continuas. De esta forma, para las variables independientes, vamos a considerar las variables cuantitativas y cualitativas que están codificadas numéricamente.

Vamos a trabajar con una partición del dataset, como conjunto de entrenamiento (con aproximadamente el 80% de los registros), sobre el que realizaremos el análisis discriminante y otra partición, como conjunto test (con el 20% de los registros restantes), con el que realizaremos la validación del modelo. Usaremos el dataset data_xlsx_G3_cathegorical definido a partir de data_xlsx_reduced (en los histogramas de la exploración gráfica de datos que se realiza a continuación se justifica el uso de este en vez del original).

# Partitioning the dataset: training (80%) + test (20%)
train_index<-createDataPartition(data_xlsx_G3_cathegorical$G3, p=0.80)$Resample1

train_data<-data_xlsx_G3_cathegorical[train_index,]
test_data<-data_xlsx_G3_cathegorical[-train_index,]

5.1 Exploración gráfica de los datos

En primer lugar exploramos como de bien (o mal) clasifica en las calificaciones cada una de las variables explicativas consideradas de forma independiente. Para ello, dibujamos los histogramas superpuestos. Si los histogramas se separan, la variable considerada sería un buen clasificador individual para la calificación.

# Variables list
variables <- c("sex", "age", "famsize", "Mjob", "reason", "traveltime",
               "studytime", "famsup", "paid", "activities", "romantic",
               "famrel", "freetime", "goout", "Dalc", "Walc", "health",
               "absences")

# Graphics list
graphics <- list()

for (var in variables) {
  p <- ggplot(data = train_data, aes_string(x = var, fill = "G3")) +
       geom_bar(position = "identity", alpha = 0.5)
  graphics[[var]] <- p
}

ggarrange(plotlist = graphics, nrow = 5, ncol = 4, common.legend = TRUE)

En este sentido, a primera vista no se ve ninguna variable que diferencia muy bien la calificación.

Si hubiéramos usado el dataset data_xlsx, hubiésemos visto que hay histogramas de algunas variables con una sola barra. Esto es porque dichas variables tienen el mismo valor para todos los registros (consecuencia del tratamiento de outliers), haciendo que su varianza sea cero (lo cual nos daría problemas para el Test de normalidad univariante de Shapiro-Wilks, el cual realizaremos posteriormente). En concreto, estas variables son school, address, Pstatus, failures, schoolsup, nursery, higher e internet. Por ello, como tenemos un gran número de variables, se ha decidido eliminarlas, tal y como hicimos en ACP y AF, resultando así el dataset que ya hemos definido anteriormente data_xlsx_reduced (el cual estamos utilizando).

A continuación, se explora qué pares de variables separan mejor entre calificaciones.

pairs(x = train_data[, c("sex", "age", "famsize", "Mjob", "reason", "traveltime",
                         "studytime", "famsup", "paid", "activities", "romantic",
                         "famrel", "freetime", "goout", "Dalc", "Walc", "health",
                         "absences")
                    ],
      col = c("firebrick", "green3")[data_xlsx_G3_cathegorical$G3], pch = 19)

Vemos que al disponer de un número de variables elevado, no se aprecia apenas el gráfico, por lo que no podemos deducir a primera vista alguna pareja de variables que separe adecuadamente las calificaciones. Igualmente, no parece que la haya.

5.2 Normalidad univariante y multivariante

5.2.1 Normalidad univariante

A continuación se realiza una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.

Histogramas univariantes

# Histogram representation of each variable for each species
par(mfcol = c(2, 3))
for (k in 1:18) {
  j0 <- names(train_data)[k]
  x0 <- seq(min(train_data[, k]), max(train_data[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(train_data$G3)[i]
    x <- train_data[train_data$G3 == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("G3", i0), xlab = j0)
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)
  }
}

Gráficos qqplots

# Representation of normal quantiles of each variable for each species
par(mfcol = c(2, 3))
for (k in 1:18) {
  j0 <- names(train_data)[k]
  x0 <- seq(min(train_data[, k]), max(train_data[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(train_data$G3)[i]
    x <- train_data[train_data$G3 == i0, j0]
    qqnorm(x, main = paste("G3", i0, j0), pch = 19, col = i + 1)
    qqline(x)
  }
}

Este análisis exploratorio puede darnos una idea de la posible distribución normal de las variables univariadas, pero siempre es mejor hacer los respectivos test de normalidad.

Test de normalidad univariantes (Shapiro-Wilks)

Se contrasta la hipótesis nula de que los datos siguen una distribución normal univariante. Se rechaza esta hipótesis si el p-valor dado por el test de Shapiro-Wilk es menor que 0.05. En otro caso no se rechaza el supuesto de normalidad de los datos.

data_tidy <- melt(train_data, value.name = "valor")
aggregate(valor ~ G3 + variable, data = data_tidy,
          FUN = function(x){shapiro.test(x)$p.value})

Se encuentran evidencias de falta de normalidad univariante (p-valor < 0.05). A pesar de ello, se procede como en ACP y AF, es decir, tiramos para adelante con el fin de mostrar cómo se realizaría el análisis discriminante, pero habría que tener cuidado a la hora de obtener los resultados.

5.2.2 Normalidad multivariante

Test de normalidad multivariante (Mardia, Henze-Zirkler y Royston)

El paquete MVN contiene funciones que permiten realizar los tres test que se utilizan habitualmente para contrastar la normalidad multivariante.

Esta normalidad multivariante puede verse afectada por la presencia de outliers multivariantes. En este paquete también encontramos funciones para el análisis de outliers multivariantes. Se haría de la siguiente forma:

\(outliers <- mvn(data = train\_data[, -ncol(train\_data)], mvnTest = "hz", multivariateOutlierMethod = "quan")\)

Si la ejecutamos obtenemos el error de que la matriz de correlaciones no tiene inversa ya que tiene determinante igual a cero, lo cual no es cierto como se puede ver a continuación.

cor_m <- cor(train_data[, -ncol(train_data)])
det(cor_m)
## [1] 0.2048836

De igual forma, los dos test realizados a continuación encuentran evidencias al 5% de significación de falta de normalidad multivariante.

# Royston multivariate normality test
royston_test <- mvn(data = train_data[,-ncol(train_data)], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality
# Henze-Zirkler multivariate normality test
hz_test <- mvn(data = train_data[,-ncol(train_data)], mvnTest = "hz")
hz_test$multivariateNormality

Cabe destacar que el discriminante lineal es muy sensible a la falta de normalidad multivariante y habría que leer sus resultados con cuidado, pero el cuadrático es bastante robusto frente a esta ausencia.

5.3 Homogeneidad de la varianza

Para el análisis de la homogeneidad de la varianza:

  • Cuando hay un solo predictor el test más recomendable es el test de Barttlet.

  • Cuando se emplean múltiples predictores, se tiene que contrastar que la matriz de covarianzas es constante en todos los grupos. En este caso es también recomendable comprobar la homogeneidad de la varianza para cada predictor a nivel individual. El test más recomendable es el test de Box M, que es una extensión del de Barttlet para escenarios multivariantes. Hay que tener en cuenta que es muy sensible a que los datos efectivamente se distribuyan según una normal multivariante. Por este motivo se recomienda utilizar una significación (p-value) < 0.001 para rechazar la hipótesis nula.

La hipóstesis nula a contrastar es la de igualdad de matrices de covarianzas en todos los grupos.

boxM(data = train_data[, -19], grouping = train_data[, 19])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  train_data[, -19]
## Chi-Sq (approx.) = 174.12, df = 171, p-value = 0.4192

En este caso no se rechaza la hipótesis nula ya que p-valor > 0.001 y por tanto se asume la homogeneidad de varianzas. Es importante recordar que para que esta conclusión sea fiable debe darse el supuesto de normalidad multivariante (la cual vimos que no se daba en nuestro caso, pero aún así tiramos para adelante).

5.4 Función discriminante

Aunque no se satisfagan los supuestos de normalidad multivariante, vamos a ajustar un modelo de clasificación discriminante lineal (tiramos para adelante con el fin de mostrar cómo se realizaría el análisis discriminante lineal, aunque hay que tenerlo presente ante la posiblidad de obtener resultados inesperados).

La función lda() del paquete MASS ajusta este modelo.

lda_model <- lda(formula = G3 ~ sex + age + famsize + Mjob + reason + traveltime +
                                studytime + famsup + paid + activities + romantic + 
                                famrel + freetime + goout + Dalc + Walc + health +
                                absences,
                 data = train_data)
lda_model
## Call:
## lda(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime + 
##     famsup + paid + activities + romantic + famrel + freetime + 
##     goout + Dalc + Walc + health + absences, data = train_data)
## 
## Prior probabilities of groups:
##  aprobado  suspenso 
## 0.6708861 0.3291139 
## 
## Group means:
##                sex      age   famsize     Mjob   reason traveltime studytime
## aprobado 0.4952830 16.55660 0.7216981 2.469297 1.301887   1.377249  2.073822
## suspenso 0.4615385 17.04503 0.7211538 2.721154 1.288462   1.405635  2.070545
##             famsup      paid activities  romantic   famrel freetime    goout
## aprobado 0.3962264 0.5141509  0.4811321 0.6981132 4.157258 3.349483 3.009434
## suspenso 0.3557692 0.6442308  0.5000000 0.5769231 4.006931 3.395100 3.432692
##              Dalc     Walc   health absences
## aprobado 1.336595 2.330189 3.551887 3.907037
## suspenso 1.432896 2.480769 3.740385 4.229504
## 
## Coefficients of linear discriminants:
##                     LD1
## sex        -0.158424951
## age         0.387186233
## famsize    -0.015049772
## Mjob        0.162799150
## reason     -0.048869654
## traveltime  0.069702592
## studytime  -0.882737708
## famsup     -0.711541738
## paid        0.850381465
## activities  0.003553655
## romantic   -0.529534163
## famrel     -0.543551569
## freetime   -0.065047956
## goout       0.517370976
## Dalc        0.370985126
## Walc       -0.114904637
## health      0.193706203
## absences   -0.006099256

La salida de este objeto, muestra las probabilidades a priori de cada grupo, en este caso 0.6708861 para aprobado y 0.3291139 para suspenso, las medias de cada regresor por grupo y los coeficientes del modelo de clasificación discriminante lineal.

Una vez construido el clasificador, se pueden clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict(). Por ejemplo, vamos a clasificar todas las observaciones del dataset test.

new_observations <- test_data

predict(object = lda_model, newdata = new_observations)
## $class
##  [1] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
##  [9] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [17] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [25] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [33] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [41] aprobado aprobado aprobado aprobado aprobado suspenso aprobado aprobado
## [49] aprobado aprobado aprobado aprobado aprobado suspenso aprobado aprobado
## [57] aprobado aprobado aprobado aprobado aprobado aprobado suspenso aprobado
## [65] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## Levels: aprobado suspenso
## 
## $posterior
##      aprobado   suspenso
## 1   0.5287203 0.47127972
## 12  0.8322448 0.16775522
## 21  0.9583155 0.04168449
## 27  0.8581690 0.14183098
## 29  0.7140187 0.28598133
## 30  0.6233143 0.37668572
## 32  0.8305574 0.16944257
## 37  0.8448085 0.15519145
## 41  0.3992119 0.60078808
## 50  0.6752142 0.32478580
## 51  0.7367466 0.26325341
## 55  0.9487876 0.05121244
## 56  0.8248841 0.17511593
## 57  0.9180837 0.08191632
## 59  0.7584683 0.24153170
## 77  0.8855864 0.11441364
## 86  0.7341848 0.26581517
## 92  0.8983315 0.10166847
## 96  0.8713945 0.12860552
## 100 0.6920912 0.30790881
## 102 0.6052546 0.39474545
## 104 0.7184537 0.28154634
## 106 0.7932028 0.20679716
## 107 0.9241282 0.07587181
## 113 0.6925321 0.30746786
## 114 0.9048116 0.09518845
## 115 0.8839749 0.11602513
## 117 0.7371635 0.26283651
## 124 0.5809519 0.41904815
## 131 0.6136074 0.38639265
## 137 0.5518007 0.44819929
## 139 0.6277175 0.37228252
## 141 0.8569860 0.14301402
## 143 0.8899853 0.11001465
## 152 0.5356066 0.46439336
## 156 0.9037751 0.09622488
## 160 0.5052616 0.49473838
## 166 0.8369355 0.16306454
## 176 0.7931131 0.20688693
## 185 0.8789448 0.12105515
## 197 0.7963019 0.20369812
## 199 0.5084081 0.49159186
## 204 0.7727307 0.22726929
## 210 0.7688319 0.23116814
## 211 0.6156840 0.38431599
## 221 0.2935173 0.70648271
## 223 0.9015213 0.09847869
## 237 0.8348371 0.16516295
## 257 0.7818562 0.21814376
## 260 0.8780754 0.12192462
## 263 0.8102649 0.18973507
## 265 0.5543106 0.44568938
## 275 0.6878363 0.31216375
## 276 0.4653202 0.53467976
## 280 0.8891487 0.11085128
## 286 0.8087309 0.19126909
## 288 0.7580155 0.24198453
## 289 0.8025402 0.19745977
## 294 0.7617490 0.23825104
## 295 0.8384137 0.16158628
## 296 0.7966999 0.20330005
## 297 0.5492399 0.45076007
## 298 0.3752390 0.62476097
## 300 0.8092060 0.19079400
## 302 0.8499428 0.15005725
## 310 0.4715284 0.52847161
## 316 0.5394510 0.46054901
## 329 0.8247123 0.17528768
## 339 0.7196361 0.28036391
## 340 0.6953299 0.30467010
## 343 0.5568235 0.44317647
## 344 0.5309468 0.46905321
## 354 0.1935134 0.80648664
## 358 0.5618232 0.43817679
## 365 0.6670069 0.33299308
## 367 0.7672386 0.23276139
## 379 0.7346458 0.26535417
## 388 0.7515503 0.24844974
## 389 0.7183630 0.28163704
## 
## $x
##              LD1
## 1    0.886774586
## 12  -0.982390749
## 21  -2.910417203
## 27  -1.232026790
## 29  -0.119059037
## 30   0.398135389
## 32  -0.967255462
## 37  -1.099108605
## 41   1.545311460
## 50   0.111181644
## 51  -0.262576290
## 55  -2.639027832
## 56  -0.917228217
## 57  -2.007077792
## 59  -0.407387406
## 77  -1.541663112
## 86  -0.246020608
## 92  -1.708124085
## 96  -1.374331598
## 100  0.013046932
## 102  0.493983864
## 104 -0.146496000
## 106 -0.658904361
## 107 -2.111706963
## 113  0.010444219
## 114 -1.799967316
## 115 -1.521787377
## 117 -0.265280286
## 124  0.620629736
## 131  0.449860110
## 137  0.769916606
## 139  0.374500688
## 141 -1.219848124
## 143 -1.597189067
## 152  0.851996350
## 156 -1.784910111
## 160  1.004914006
## 166 -1.025115103
## 176 -0.658216330
## 185 -1.461251985
## 197 -0.682792002
## 199  0.989086100
## 204 -0.507338589
## 210 -0.479591758
## 211  0.438836347
## 221  2.135769816
## 223 -1.752660679
## 237 -1.005881887
## 257 -0.573627127
## 260 -1.451009218
## 263 -0.793931666
## 265  0.757149679
## 275  0.038056697
## 276  1.206073713
## 280 -1.586481128
## 286 -0.781424235
## 288 -0.404281403
## 289 -0.731712057
## 294 -0.430009130
## 295 -1.038784224
## 296 -0.685879859
## 297  0.782928478
## 298  1.672372407
## 300 -0.785289623
## 302 -1.149026703
## 310  1.174725459
## 316  0.832552201
## 329 -0.915733887
## 339 -0.153855209
## 340 -0.006118322
## 343  0.744353355
## 344  0.875536714
## 354  2.826017934
## 358  0.718849065
## 365  0.157935977
## 367 -0.468347528
## 379 -0.248992315
## 388 -0.360359713
## 389 -0.145932280

Por ejemplo, según la función discriminante, la probabilidad a posteriori de que el segundo registro del conjunto test apruebe es del 89.6% mientras la de que suspenda es tan solo del 10.4%. Por tanto este dato sería clasificado en la calificación aprobado. Se razona del mismo modo con el resto de elementos del conjunto test.

5.5 Validación del modelo

La función confusionmatrix() del paquete biotools realiza una validación cruzada del modelo de clasificación.

pred <- predict(object = lda_model, newdata = test_data)
confusionmatrix(test_data$G3, pred$class)
##          new aprobado new suspenso
## aprobado           50            3
## suspenso           23            3
# Classification error percentage
test_error <- mean(test_data$G3 != pred$class) * 100
paste("test_error=", test_error, "%")
## [1] "test_error= 32.9113924050633 %"

5.6 Visualización de las clasificaciones

Desde el punto de vista geométrico, el análisis discriminante lineal separa el espacio mediante una recta. En este sentido, la función partimat() del paquete klaR permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.

partimat(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
              famsup + paid + activities + romantic + famrel + freetime + 
              goout + Dalc + Walc + health + absences,
         data = test_data, method = "lda", prec = 200,
         image.colors = c("green", "orange"),
         col.mean = "yellow",nplots.vert =1, nplots.hor=3)

6 Análisis Discriminante Cuadrático (ADC)

Al igual que para el Análisis Discriminante Lineal, para realizar el Análisis Discriminante Cuadrático se empieza con la exploración gráfica de los datos y las comprobaciones sobre normalidad univariante y multivariante y homogeneidad de las varianzas, que ya se han realizado previamente.

Vimos que no se verifica el supuesto de normalidad multivariante. Aun así, se procede a ajustar un modelo discriminante cuadrático porque es robusto frente a la falta de este supuesto, aunque hay que tenerlo presente ante la posiblidad de obtener resultados inesperados.

6.1 Función discriminante

La función qda() del paquete MASS realiza la clasificación.

qda_model <- qda(formula = G3 ~ sex + age + famsize + Mjob + reason + traveltime +
                                studytime + famsup + paid + activities + romantic + 
                                famrel + freetime + goout + Dalc + Walc + health +
                                absences,
                 data = train_data)
qda_model
## Call:
## qda(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime + 
##     famsup + paid + activities + romantic + famrel + freetime + 
##     goout + Dalc + Walc + health + absences, data = train_data)
## 
## Prior probabilities of groups:
##  aprobado  suspenso 
## 0.6708861 0.3291139 
## 
## Group means:
##                sex      age   famsize     Mjob   reason traveltime studytime
## aprobado 0.4952830 16.55660 0.7216981 2.469297 1.301887   1.377249  2.073822
## suspenso 0.4615385 17.04503 0.7211538 2.721154 1.288462   1.405635  2.070545
##             famsup      paid activities  romantic   famrel freetime    goout
## aprobado 0.3962264 0.5141509  0.4811321 0.6981132 4.157258 3.349483 3.009434
## suspenso 0.3557692 0.6442308  0.5000000 0.5769231 4.006931 3.395100 3.432692
##              Dalc     Walc   health absences
## aprobado 1.336595 2.330189 3.551887 3.907037
## suspenso 1.432896 2.480769 3.740385 4.229504

La salida de este objeto, muestra las probabilidades a priori de cada grupo, en este caso 0.6708861 para aprobado y 0.3291139 para suspenso, las medias de cada regresor por grupo.

Una vez construido el clasificador, se pueden clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict(). Por ejemplo, vamos a clasificar todas las observaciones del dataset test.

new_observations <- test_data

predict(object = qda_model, newdata = new_observations)
## $class
##  [1] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
##  [9] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [17] aprobado suspenso aprobado aprobado aprobado suspenso aprobado aprobado
## [25] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
## [33] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [41] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [49] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [57] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [65] aprobado aprobado aprobado aprobado aprobado suspenso aprobado suspenso
## [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## Levels: aprobado suspenso
## 
## $posterior
##      aprobado     suspenso
## 1   0.9420895 5.791050e-02
## 12  0.9132794 8.672064e-02
## 21  0.9983882 1.611764e-03
## 27  0.9996942 3.058088e-04
## 29  0.9519198 4.808024e-02
## 30  0.8740865 1.259135e-01
## 32  0.9725168 2.748321e-02
## 37  0.9933374 6.662618e-03
## 41  0.1599467 8.400533e-01
## 50  0.3660455 6.339545e-01
## 51  0.9921764 7.823598e-03
## 55  0.9978306 2.169449e-03
## 56  0.9628626 3.713736e-02
## 57  0.9866293 1.337075e-02
## 59  0.9659756 3.402435e-02
## 77  0.7117161 2.882839e-01
## 86  0.6263690 3.736310e-01
## 92  0.4917608 5.082392e-01
## 96  0.7509006 2.490994e-01
## 100 0.8578510 1.421490e-01
## 102 0.6992017 3.007983e-01
## 104 0.2780915 7.219085e-01
## 106 0.9881970 1.180300e-02
## 107 0.9652249 3.477514e-02
## 113 0.5679399 4.320601e-01
## 114 0.9961207 3.879346e-03
## 115 0.9992368 7.632085e-04
## 117 0.9881334 1.186660e-02
## 124 0.1060614 8.939386e-01
## 131 0.1369429 8.630571e-01
## 137 0.9139531 8.604690e-02
## 139 0.5283120 4.716880e-01
## 141 0.7562291 2.437709e-01
## 143 0.9327727 6.722733e-02
## 152 0.6146317 3.853683e-01
## 156 0.9999709 2.906058e-05
## 160 0.4715576 5.284424e-01
## 166 0.8308611 1.691389e-01
## 176 0.6283201 3.716799e-01
## 185 0.3972366 6.027634e-01
## 197 0.9960218 3.978206e-03
## 199 0.8719620 1.280380e-01
## 204 0.9873704 1.262960e-02
## 210 0.7181323 2.818677e-01
## 211 0.5816051 4.183949e-01
## 221 0.5186015 4.813985e-01
## 223 0.9791096 2.089043e-02
## 237 0.8131715 1.868285e-01
## 257 0.7129042 2.870958e-01
## 260 0.9856287 1.437128e-02
## 263 0.9499050 5.009498e-02
## 265 0.7846896 2.153104e-01
## 275 0.7546396 2.453604e-01
## 276 0.9453512 5.464877e-02
## 280 0.9917348 8.265192e-03
## 286 0.9869410 1.305902e-02
## 288 0.5607798 4.392202e-01
## 289 0.9415090 5.849101e-02
## 294 0.6537919 3.462081e-01
## 295 0.8604681 1.395319e-01
## 296 0.6809333 3.190667e-01
## 297 0.9434945 5.650552e-02
## 298 0.7037053 2.962947e-01
## 300 0.8192494 1.807506e-01
## 302 0.8940954 1.059046e-01
## 310 0.8961285 1.038715e-01
## 316 0.5901816 4.098184e-01
## 329 0.9292400 7.076001e-02
## 339 0.9845709 1.542907e-02
## 340 0.4540543 5.459457e-01
## 343 0.8105703 1.894297e-01
## 344 0.4402601 5.597399e-01
## 354 0.1499105 8.500895e-01
## 358 0.8961398 1.038602e-01
## 365 0.8286849 1.713151e-01
## 367 0.9680180 3.198197e-02
## 379 0.6149508 3.850492e-01
## 388 0.9583831 4.161691e-02
## 389 0.9928296 7.170400e-03

Por ejemplo, según la función discriminante, la probabilidad a posteriori de que el segundo registro del conjunto test apruebe es del 99.9% mientras la de que suspenda es tan solo del 0.1%. Por tanto este dato sería clasificado en la calificación aprobado. Se razona del mismo modo con el resto de elementos del conjunto test.

6.2 Validación del modelo

La función confusionmatrix() del paquete biotools realiza una validación cruzada del modelo de clasificación.

pred <- predict(object = qda_model, newdata = test_data)
confusionmatrix(test_data$G3, pred$class)
##          new aprobado new suspenso
## aprobado           47            6
## suspenso           21            5
# Classification error percentage
test_error <- mean(test_data$G3 != pred$class) * 100
paste("test_error=", test_error, "%")
## [1] "test_error= 34.1772151898734 %"

6.3 Visualización de las clasificaciones

La función partimat() del paquete klaR permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.

partimat(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
              famsup + paid + activities + romantic + famrel + freetime + 
              goout + Dalc + Walc + health + absences,
         data = test_data, method = "qda", prec = 200,
         image.colors = c("darkgoldenrod1", "skyblue2"),
         col.mean = "firebrick", nplots.vert =1, nplots.hor=3)

LDA vs QDA

El clasificador más adecuado depende de las implicaciones que tenga asumir que todos los grupos comparten una matriz de covarianza común ya que este supuesto puede producir un sesgo en las clasificaciones o producir varianzas altas.

  • LDA produce límites de decisión lineales, lo que se traduce en menor flexibilidad y por lo tanto menor problema de varianza.
  • QDA produce límites cuadráticos y por lo tanto curvos, lo que aporta mayor flexibilidad permitiendo ajustarse mejor a los datos, menor sesgo pero mayor riesgo de varianza.
  • En términos generales, LDA tiende a conseguir mejores clasificaciones que QDA cuando hay pocas observaciones con las que entrenar al modelo, escenario en el que evitar la varianza es crucial.
  • Si se dispone de una gran cantidad de observaciones de entrenamiento (como es nuestro caso) o si no es asumible que existe una matriz de covarianza común entre clases, QDA es más adecuado.

7 Análisis cluster (AC)

Vamos a realizar un Análisis Clúster que confirme que la agrupación de la variable respuesta (variable G3 categórica) utilizada en los modelos de clasificación lineal y cuadrático es adecuada.

7.1 Algunas distancias para Análisis Cluster

Para clasificar las observaciones en grupos es necesario elegir medidas de similaridad, o de distancia (disimilaridad), adecuadas que proporcionen información de como de parecidas son dos observaciones cualesquiera. De hecho esta elección influye en el tamaño y forma de los clusters. Esta elección es un paso fundamental en clustering.

Algunas medidas de distancia clásicas frecuentemente utilizadas son la distancia Euclídea o la distancia Manhattan.

Otras medidas de disimilaridad ampliamente utilizadas, por ejemplo, en el análisis de datos de expresión génica son las distancias basadas en correlaciones. Estas distancias se obtienen restando la correspondiente medida de correlación al valor 1. Entre estas medidas de distancia destacan: la distancia basada en la correlación de Pearson, en la correlación de Spearman, correlación de Kendall, etc.

Como se ha dicho anteriormente la elección de la distancia es muy importante. Casi todo el software habitual para Análisis Cluster utiliza la distancia euclídea, aunque dependiendo del tipo de datos y de las preguntas de investigación planteadas puede interesar otra medida de disimilaridad o distancia.

En R es fácil calcular y visualizar la matriz de distancias entre observaciones con las funciones get_dist() y fviz_dist(), respectivamente, incluidas en el paquete factoextra. De hecho, aunque por defecto, get_dist() calcula la distancia euclídea también admite como parámetro todas las distancias mencionadas anteriormente.

La siguiente matriz de distancias muestra en marrón aquellos estados que presentan grandes disimilitudes (distancias), frente a aquellos que parecen más cercanos en amarillo. Se utiliza el color blanco para referirse a aquellos estados con distancias no tan extremas como para ser consideradas como bajas o altas.

distance<- get_dist(data_xlsx_reduced)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))

7.2 Clustering jerárquico: método de Ward

El agrupamiento jerárquico está interesado en encontrar una jerarquía basada en la cercanía o semejanza de los datos según la distancia considerada. En el caso aglomerativo, se parte de un grupo con las observaciones más cercanas. A continuación se calculan los siguientes pares más cercanos y de manera ascendente se van generando grupos. Esta construcción se puede observar de forma visual con la construcción de un dendrograma.

A continuación se ilustrará cómo los grupos están definidos por la cantidad de líneas verticales del dendrograma, y la selección del número de grupos óptimo se podrá estimar desde este mismo gráfico.

dendrogram <- hclust(dist(data_xlsx_reduced, method = 'euclidean'), method = 'ward.D')
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) + 
             labs(title = "Dendrograma")

En el eje horizontal del dendrograma tenemos cada uno de los datos que componen el conjunto de entrada, mientras que en el eje vertical se representa la distancia euclídea que existe entre cada grupo a medida que éstos se van jerarquizando.

Cada línea vertical del diagrama representa un agrupamiento. Conforme se va subiendo en el dendograma termina con un solo gran grupo determinado por la línea horizontal superior. De modo que, al ir descendiendo en la jerarquía, se observa que de un solo grupo pasamos a 2, luego a 3, luego a 4, y así sucesivamente.

Una forma de determinar el número K de grupos adecuado es cortando el dendrograma a aquella altura del diagrama que mejor representa los datos de entrada.

7.3 Clustering no jerárquico: algoritmo K-means

R implementa el algoritmo K-means con la función del mismo nombre. Esta función recibe como parámetros de entrada los datos y el número de agrupamientos a realizar (parámetro centers). Para abordar el problema de la elección de los puntos semilla iniciales incorpora el parámetro nstart que prueba múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, si nstart = 25, generará 25 configuraciones iniciales. El uso de este parámetro es recomendable.

Para este primer ejemplo la función kmeans() construye dos clusters.

k2 <- kmeans(data_xlsx_reduced, centers = 2, nstart = 25)

# Displaying all the fields of the object k2
str(k2)
## List of 9
##  $ cluster     : int [1:395] 1 1 1 1 1 1 1 1 1 1 ...
##  $ centers     : num [1:2, 1:19] 0.484 0.386 16.646 16.977 0.698 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "1" "2"
##   .. ..$ : chr [1:19] "sex" "age" "famsize" "Mjob" ...
##  $ totss       : num 19315
##  $ withinss    : num [1:2] 13240 750
##  $ tot.withinss: num 13990
##  $ betweenss   : num 5325
##  $ size        : int [1:2] 351 44
##  $ iter        : int 1
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

La salida que proporciona la función kmeans() es una lista de información, de la que destacan las siguientes:

  • cluster: es un vector de enteros, de 1 a K (K = 2 en este caso), que indica el cluster en el que ha sido ubicado cada observación.
  • centers: una matriz con los sucesivos centros de los clusters.
  • totss: la suma total de cuadrados.
  • withinss: vector de suma de cuadrados dentro de cada cluster (un componente por cluster).
  • tot.withinss: suma total de cuadrados de los cluster, i.e. sum(withinss).
  • betweens: suma de cuadrados entre grupos, i.e. totss-tot.withinss.
  • size: el número de observaciones en cada cluster.

Al mostrar la variable k2 se ve como las agrupaciones dan como resultado 2 tamaños de agrupación de 44 y 351 También se ven los centros de cada grupo (medias) en todas las variables. Y por último la asignación de grupo para cada observación.

k2
## K-means clustering with 2 clusters of sizes 351, 44
## 
## Cluster means:
##         sex      age   famsize     Mjob   reason traveltime studytime    famsup
## 1 0.4843305 16.64582 0.6980057 2.458652 1.270655   1.379977  2.072080 0.3874644
## 2 0.3863636 16.97727 0.8181818 2.818182 1.295455   1.517970  2.076143 0.3863636
##        paid activities  romantic   famrel freetime    goout     Dalc     Walc
## 1 0.5099715  0.4985755 0.6894587 4.109019 3.353057 3.079772 1.333515 2.313390
## 2 0.7954545  0.4318182 0.4772727 4.054816 3.311291 3.340909 1.363636 2.113636
##     health  absences       G3
## 1 3.535613 4.4491379 11.62678
## 2 3.704545 0.2860918  0.75000
## 
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 1 2 1 1 1 2 1 2 1
## [149] 2 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
## [334] 2 2 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 13240.3490   750.1396
##  (between_SS / total_SS =  27.6 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Una forma visual de resumir los resultados de forma elegante y con una interpretación directa es mediante el uso de la función fviz_cluster().

fviz_cluster(k2, data=data_xlsx_reduced)

Observación: Si hay más de dos dimensiones (variables), esta función realizará, en primer lugar, un análisis de componentes principales (ACP) y dibujará los puntos de acuerdo a las dos primeras componentes principales obtenidas (las que explican la mayor parte de la varianza). Es por esto que en el gráfico anterior, Dim1 y Dim2 se refieren a estas dos componentes principales.

Para usar el algoritmo K-means, el número K de clusters debe ser fijado de antemano. Es por esto que es recomendable ejecutar el mismo proceso con otros valores de K (en este ejemplo para K = 3, 4 y 5) para comparar y examinar las diferencias entre los resultados.

set.seed(123)

k3 <- kmeans(data_xlsx_reduced, centers = 3, nstart = 25)
k4 <- kmeans(data_xlsx_reduced, centers = 4, nstart = 25)
k5 <- kmeans(data_xlsx_reduced, centers = 5, nstart = 25)

# Plots to compare
p24 <- fviz_cluster(k2, geom = "point", data = data_xlsx_reduced) + ggtitle("k = 2")
p25 <- fviz_cluster(k3, geom = "point",  data = data_xlsx_reduced) + ggtitle("k = 3")
p26 <- fviz_cluster(k4, geom = "point",  data = data_xlsx_reduced) + ggtitle("k = 4")
p27 <- fviz_cluster(k5, geom = "point",  data = data_xlsx_reduced) + ggtitle("k = 5")


grid.arrange(p24, p25, p26, p27, nrow = 2)

Aunque esta visualización nos permite deducir donde ocurren las verdaderas diferencias (o no ocurren, ya que vemos que los clusters se superponen), no nos dice cuál es el número óptimo de clusters.

7.3.1 Determinación del número óptimo de clusters

Tal y como se ha indicado anteriormente, cuando se aplica un método no jerárquico como K-means para realizar un análisis cluster, el investigador debe informar a priori del número de clusters deseado. En este sentido, este investigador estará interesado en proporcionar de partida un número óptimo de grupos a formar.

A continuación se presentan tres de los métodos más utilizados para determinar este número óptimo de grupos: método de Elbow, método de Silhouette y el stadístico Gap.

7.3.1.1 Método de Elbow

Teniendo en mente que la idea tras una división en K clusters, es obtener estas agrupaciones de modo que la varianza total intra-grupos sea mínima (total within-cluster variation o total within-cluster sum of square), se puede utilizar el siguiente algoritmo para identificar el número óptimo de clusters:

  • Ejecutar un algoritmo de clustering (como K-means) para diferentes valores del K (por ejemplo K = 1,…,10).
  • Para cada K se calcula la variación total intra-cluster (total within-cluster sum of square, que aquí denotamos por wss).
  • Se dibuja la curva de wss de acuerdo al número de clusters K. La localización en esta curva de un “codo” es tomado como el indicador más apropiado del número de clusters.

Aunque podemos programar este algoritmo en R, el método de Elbow está implementado en la función fviz_nbclust().

set.seed(123)
fviz_nbclust(data_xlsx_reduced, kmeans, method = "wss")

7.3.1.2 Método de Silhouette

Este enfoque, silueta promedio, mide la calidad de un cluster. Es decir, determina como de adecuado es un objeto dentro de su cluster. El número óptimo de clusters según este enfoque es, de entre un rango de valores posibles para K, aquel que maximiza la silueta promedio.

Como antes, este algoritmo se puede programar en R, pero la función fviz_nbclust() también lo incluye.

set.seed(123)
fviz_nbclust(data_xlsx_reduced, kmeans, method = "silhouette")

7.3.1.3 Método estadístico de brecha (GAP)

Este método compara la variación total intracluster para diferentes valores de K. Utiliza simulación Montecarlo en su algoritmo.

La función clusGap() proporciona el estadístico GAP y su error estándar para una salida. Con la función fviz_gap_stat() se obtiene una representación gráfica que sugiere un número de clusters.

set.seed(123)
gap_stat <- clusGap(data_xlsx_reduced, FUN = kmeans, nstart = 25, K.max = 10, B = 50)
fviz_gap_stat(gap_stat)

7.3.2 Análisis de los resultados

Tras el análisis pormenorizado del número óptimo de clusters, realizado en la sección anterior, y las pruebas anteriores con diferentes números de clusters, concluímos que no parece adecuado realizar un Análisis Cluster. Esto puede deberse quizá a la independencia de las variables que observamos a la hora de realizar ACP o AF o la forma en la que se han codificado las variables categóricas (podría probarse a hacer este análisis cluster usando medidas de distancia para variables categóricas, como puede ser la distancia de Manhattan). Se observa que a mayor número de clusters, peor (ya que los clusters se superponen)

Aun así, si tuviesemos que realizarlo, K = 2 parece ser el número de grupos más adecuado. Esto confirma que la agrupación de la variable respuesta utilizada en los modelos de clasificación lineal y cuadrático es adecuada. Por lo tanto, el categorizar la calificación de los alumnos como aprobado o suspenso, fue una buena opción tomada a la hora de realizar los modelos de clasificación (por ejemplo, podríamos haber pensado también en categorizar la calificación de los alumnos como suspenso, aprobado, notable y sobresaliente).

set.seed(123)
final <- kmeans(data_xlsx_reduced, centers = 2, nstart = 25)
print(final)
## K-means clustering with 2 clusters of sizes 351, 44
## 
## Cluster means:
##         sex      age   famsize     Mjob   reason traveltime studytime    famsup
## 1 0.4843305 16.64582 0.6980057 2.458652 1.270655   1.379977  2.072080 0.3874644
## 2 0.3863636 16.97727 0.8181818 2.818182 1.295455   1.517970  2.076143 0.3863636
##        paid activities  romantic   famrel freetime    goout     Dalc     Walc
## 1 0.5099715  0.4985755 0.6894587 4.109019 3.353057 3.079772 1.333515 2.313390
## 2 0.7954545  0.4318182 0.4772727 4.054816 3.311291 3.340909 1.363636 2.113636
##     health  absences       G3
## 1 3.535613 4.4491379 11.62678
## 2 3.704545 0.2860918  0.75000
## 
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 1 2 1 1 1 2 1 2 1
## [149] 2 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
## [334] 2 2 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 13240.3490   750.1396
##  (between_SS / total_SS =  27.6 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# Final output
fviz_cluster(final, data = data_xlsx_reduced)